SUBROUTINE TYPE67 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C************************************************************************ C* SUBROUTINE: ACDX 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* outlet air temperature and humidity, total C* and sensible cooling capacities and power C* consumption. Performance does NOT include C* the effects of indoor fan. 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) TAirLvg Leaving air dry bulb temperature(C) 14.9211 C* OUT(2) WAirLvg Leaving air humidity ratio(-) .00876644 C* OUT(3) QTot Total heat transfer rate(W) 67198.0 C* OUT(4) QSen Sensible heat transfer rate(W) 51494.0 C* OUT(5) Power Power consumption of compressor section(W) 26945.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*********************************************************************** 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: January 1, 1992 C C INCLUDE FILES: acdx.inc C prop.inc C SUBROUTINES CALLED: BYPASS C DXDOE C TDB_TWB C FUNCTIONS CALLED: DRYBULB C ENTHALPY 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 P(MRated) Airflow at rating (kg/s) C P(UACp) Bypass parameter: BF = EXP(-UACp/Mrated) (kg/s) C hDelta Enthalpy difference (J/kg) C wDelta Humidity ratio difference (-) C mPerWRated Airflow per unit capacity at rating (kg/s W) C tAmbRated Outdoor ambient dry bulb temperature at rating (C) C twbEntRated Entering air wet bulb temperature at rating (C) C hEntRated Entering air enthalpy at rating (J/kg) C tEntRated Entering air dry bulb temperature at rating (C) C wEntRated Entering air humidity ratio at rating (-) C hLvgRated Leaving air enthalpy at rating (J/kg) C tLvgRated Leaving air dry bulb temperature at rating (C) C wLvgRated Leaving air humidity ratio at rating (-) C wAdpRated Apparatus dew point humidity ratio at rating (-) C bfRated Bypass factor at rating (-) C shr Sensible heat ratio (-) C cop Coefficient of performance (-) C wDry Humidity ratio at which SHR=1 (-) C hAirLvg Leaving air enthalpy (J/kg) 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 REAL mPerWRated,PROP(16),MRATED,UACP,PAR INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd,ERRSTAT DIMENSION PAR(3),XIN(4),OUT(5) CHARACTER*3 YCHECK(4),OCHECK(5) 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 IOPT = -1 NI = 4 !CORRECT NUMBER OF INPUTS NP = 3 !CORRECT NUMBER OF PARAMETERS ND = 0 !CORRECT NUMBER OF DERIVATIVES DATA itmax/20/ DATA YCHECK/'MF2','TE1','DM1','TE1'/ DATA OCHECK/'TE1','DM1','PW2','PW2','PW2'/ MAIR = XIN(1) TAIRENT = XIN(2) WAIRENT = XIN(3) TAIRAMB = XIN(4) CAPRATED = PAR(1) COPRATED = PAR(2) SHRRATED = PAR(3) C2*** ARI rating conditions in SI units corresponding to C2*** tAmbRated = 95 F, tEntRated = 80 F, twbEntRated = 67 F, and C2*** rated flow = 450 cfm/ton standard air DATA mPerWRated /0.00007255/ DATA tEntRated /26.6666667/ DATA twbEntRated /19.4444444/ ErrStat = 0 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 C2********************************************************************** C2 The code between these bars of asterisks is used to set internal C2 parameters and is independent of component input values. In an C2 hourly simulation, this block of code may be skipped after the C2 first call. C1*** Calculate psychrometric properties at rating point CALL TDB_TWB (Prop,tEntRated,wEntRated,RH,hEntRated,twbEntRated, & TDP,RhoD,RhoM,ErrStat) C1*** Calculate exit conditions at ARI rating point MRated = CapRated*mPerWRated hDelta = CapRated/MRated wDelta = hDelta*(1.-ShrRated)/Prop(Hfg) hLvgRated = hEntRated - hDelta wLvgRated = wEntRated - wDelta tLvgRated = DRYBULB(Prop,hLvgRated,wLvgRated) C1*** Calculate bypass factor and determine UA/cp using BF=exp(-UA/mcp) C2*** The relationship BF=exp(-UA/mcp) models the cooling coil as an C2*** enthalpy exchanger with Cmin/Cmax = 0. CALL BYPASS(Prop,tEntRated,wEntRated,tLvgRated, & wLvgRated,tAdpRated,wAdpRated, & bfRated,ErrStat) UACp = -log(bfRated)*MRated C2********************************************************************** C1*** Calculate capacity, efficiency and sensible heat ratio for given C1*** entering and outdoor air conditions using DOE 2.1D correlations CALL DXDOE(Prop,MRATED,CAPRATED,COPRATED,UACp,MAir,TAirEnt, & WAirEnt,TAirAmb,QTot,cop,shr,ErrStat) C1*** Check for dry evaporator coil IF (shr .GT. 1.) THEN C1*** Dry evaporator - Find entering humidity ratio such that SHR = 1 C2*** The correction factors in DOE 2.1 are a function of entering C2 wet bulb temperature and not entering dry bulb temperature, as C2 is appropriate for a wet coil. Since the performance of a C2 dry coil is not a function of humidity ratio, its performance C2 at a given entering air dry bulb temperature should be C2 approximately constant at all humidity ratios below the wDry, C2 defined as the humidity ratio at which SHR = 1. Since the DOE C2 correlations should apply at the wDry limit, the performance of C2 a dry coil is determined by iteratively solving for wDry using C2*** the DOE correlations. C1*** Initialize iteration parameters wDry = WAirEnt DO 100 iter=1,itmax error = shr - 1. wDry = XITERATE(wDry,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, exit loop IF (icvg .EQ. 1) GO TO 110 C1*** Calculate capacity, efficiency and sensible heat ratio for given C1*** entering and outdoor air conditions using DOE 2.1D correlations CALL DXDOE(Prop,MRATED,CAPRATED,COPRATED,UACP,MAir, & TAirEnt,wDry,TAirAmb,QTot,cop,shr,ErrStat) 100 CONTINUE C1*** Dry-out humidity ratio has not converged after maximum iterations. C1*** Print error message, set return error flag, and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN SUBOUTINE ACDX ***'/ & 1X,' Dry-out humidity ratio has not ' & 'converged after ',I2,' iterations'/) ErrStat = 1 ENDIF 110 CONTINUE C1*** Calculate outlet conditions, capacities, and power hDelta = QTot/MAir wDelta = hDelta*(1.-shr)/(Prop(Hfg)+Prop(CpVap)*TAirEnt) hAirLvg = ENTHALPY(Prop,TAirEnt,WAirEnt) - hDelta WAirLvg = WAirEnt - wDelta TAirLvg = DRYBULB(Prop,hAirLvg,WAirLvg) QSen = shr*QTot Power = QTot/cop 999 CONTINUE OUT(1) = TAIRLVG OUT(2) = WAIRLVG OUT(3) = QTOT OUT(4) = QSEN OUT(5) = POWER RETURN 1 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) SUBROUTINE BYPASS(Prop,TEnt,WEnt,TLvg,WLvg, & TAdp,WAdp,BF,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C************************************************************************ C* SUBROUTINE: BYPASS C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate apparatus dew point and bypass C* factor given entering and leaving moist C* air conditions of cooling coil. C*********************************************************************** C* INPUT VARIABLES C* TEnt Entering air temperature (C) C* WEnt Entering air humidity ratio (-) C* TLvg Leaving air temperature (C) C* WLvg Leaving air humidity ratio (-) C* C* OUTPUT VARIABLES C* TAdp Apparatus dewpoint temperature (C) C* WAdp Apparatus dewpoint humidity ratio (-) C* BF Bypass factor (-) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C*********************************************************************** C MAJOR RESTRICTIONS: None 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: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: DEWPOINT C ENTHALPY C HUMRATIO C SATPRESS C XITERATE C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C hEnt Entering air enthalpy C hLvg Leaving air enthalpy C hAdp Air enthalpy at apparatus dew point C slope Ratio temperature difference to humidity difference C between entering and leaving air states C tAdpEst Estimate of TAdp from slope 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) INTEGER ErrStat DATA itmax/20/ C1*** Iterate to determine apparatus dewpoint at which the ADP C1*** equals the temperature calculated by extending the line between C1*** entering and leaving conditions to the saturation curve C1*** Calculate "slope" of temperature vs. humidity ratio between C1*** entering and leaving states slope = (TEnt-TLvg)/(WEnt-WLvg) C1*** Initialize iteration parameters TAdp = DEWPOINT(Prop,WLvg) DO 100 iter=1,itmax C1*** Calculate apparatus dewpoint and compare with predicted value C1*** using entering conditions and slope WAdp = HUMRATIO(Prop(Patm),SATPRESS(Prop,TAdp)) TAdpEst = TEnt - slope*(WEnt-WAdp) error = TAdp-TAdpEst TAdp = XITERATE(TAdp,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, exit loop IF (icvg .EQ. 1) GO TO 110 100 CONTINUE C1*** Apparatus dewpoint has not converged after maximum iterations. C1*** Print error message, set return error flag, and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN SUBOUTINE BYPASS ***'/ & 1X,' Apparatus dewpoint has not ' & 'converged after ',I2,' iterations'/) ErrStat = 1 110 CONTINUE C1*** Calculate bypass factor from enthalpies hLvg = ENTHALPY(Prop,TLvg,WLvg) hEnt = ENTHALPY(Prop,TEnt,WEnt) hAdp = ENTHALPY(Prop,TAdp,WAdp) BF = (hLvg-hAdp)/(hEnt-hAdp) 999 RETURN END SUBROUTINE DXDOE (Prop,MRATED,CAPRATED,COPRATED,UACP, & MAir,TAirEnt,WAirEnt,TAirAmb, & QTot,Cop,Shr,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations 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 C* MAir Dry air mass flow rate (kg/s) C* TAirEnt Entering air dry bulb temperature (C) C* WAirEnt Entering air humidity ratio (-) C* TAirAmb Outdoor ambient dry bulb temperature (C) C* C* OUTPUT VARIABLES C* QTot Total heat transfer rate (W) C* Cop Coefficient of Performance (-) C* Shr Sensible heat ratio (-) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C* C* PARAMETERS C* CapRated Total capacity at ARI rating conditions (W) C* CopRated COP at ARI rating conditions (-) C* ShrRated Sensible heat ratio at ARI rating conditions (-) C* UACp Coefficient in bypass factor relationship (kg/s) C* Bypass = EXP(-UACp/MAir) C* MRated Airflow at ARI rating point (kg/s) 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************************************************************************ INTEGER ErrStat REAL MAir,mFract,MRATED,UACP,CAPRATED,COPRATED 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) C* DOE 2.1 Correction factor statement functions CAP_T_COR(twbF,tAmbF) = 0.6003404 + 0.00228726*twbF & - 0.0000128*twbF*twbf + 0.00138975*tAmbF & - 0.0000806*tambF*tAmbF + 0.00014125*twbF*tAmbF CAP_FLOW_COR(mFract) = 0.8 + 0.2*mFract EIR_T_COR(twbF,tAmbF) = -0.9617787 + 0.04817751*twbF & - 0.0002311*twbF*twbF + 0.00324392*tAmbF & + 0.00014876*tAmbF*tAmbF - 0.00029545*twbF*tAmbF EIR_FLOW_COR(mFract) = 1.156 - 0.1816*mFract DATA itmax/20/ ErrStat = 0 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(twbF,tAmbF)*CAP_FLOW_COR(mFract) Cop = CopRated/EIR_T_COR(twbF,tAmbF)/EIR_FLOW_COR(mFract) 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 RETURN END SUBROUTINE TDB_TWB (Prop,TDB,W,RH,H,TWB,TDP,RhoD,RhoM,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: TDB_TWB C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate psychrometric properties of C* moist air given dry bulb temperature C* and wet bulb temperature. C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* TWB Wet bulb temperature (C) C* C* OUTPUT VARIABLES C* W Humidity ratio (-) C* RH Relative humidity (-) C* H Enthalpy of air (J/kg) 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* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* CpWat Specific heat of liquid water (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) 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 RELHUM C ENTHALPY C DEWPOINT C RHODRY C RHOMOIST C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C psat Saturated water vapor pressure (Pa) C wstar Saturated humidity ratio at wet bulb temp. (-) 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 for a given wet bulb C1*** temperature. psat = SATPRESS (Prop,TWB) C1*** Calculate the humidity ratio corresponding to saturation C1*** pressure for a given wet bulb temperature. wstar = HUMRATIO (Prop(Patm),psat) C1*** Calculate the humidity ratio as a function of dry bulb C1*** temperature, wet bulb temperature and the humidity ratio C1*** at saturation pressure for the given wet bulb temperature. W = ((Prop(Hfg)+(Prop(CpVap)-Prop(CpWat))*TWB)*wstar- & Prop(CpAir)*(TDB-TWB))/ & (Prop(Hfg)+Prop(CpVap)*TDB-Prop(CpWat)*TWB) C1*** Calculate the saturation pressure for a given temperature. psat = SATPRESS (Prop,TDB) C1*** Calculate the relative humidity as a function of partial C1*** water vapor pressure and water vapor pressure at saturation. 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 dewpoint temperature as a function of 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 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 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 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