/*********************************************************************/ /* CONTENT OF PSYCHRO?.NMF LIBRARYFILES */ /* */ /* DATE: April 10, 1997 rev 1.0 */ /* Revisions: */ /* 99 MV Add PSYCHRO4 */ /* */ /* 980126 AB Add comments */ /* */ /* 971125 AB Remove side effect from ENTHAL */ /* */ /* FILES CONTAIN FOLLOWING FUNCTION AND SUBROUTINES */ /* AND THEIR HELP FILES TO GENERATE ANALYTICAL JACOBIANS. */ /* */ /* DEWPNT (W) THESE FILES ARE IN PSYCHRO1.NMF */ /* DRYBULB (H,W) */ /* ENTHAL (TDB,W) */ /* ENTHSAT (TDB) */ /* */ /* DEWPNTJ */ /* DRYBULBJ */ /* ENTHALJ */ /* ENTHSATJ */ /* */ /* HUMRAT (Patm,Pw) THESE FILES ARE IN PSYCHRO2.NMF */ /* HUMTH (TDB,H) */ /* RELHUM (Patm,Psat,HumRation) */ /* RHODRY (TDB,W) */ /* RHOMOIS (TDB,W) */ /* */ /* HUMRATJ */ /* HUMTHJ */ /* RELHUMJ */ /* RHODRYJ */ /* RHOMOISJ */ /* */ /* SATPRES (T) THESE FILES ARE IN PSYCHRO3.NMF */ /* SATTEMP (P) */ /* TAIRSAT (HSat) */ /* WETBULB (TDB,W) */ /* */ /* SATPRESJ */ /* SATTEMPJ */ /* TAIRSATJ */ /* WETBULBJ */ /* */ /* DEWPNTP (W,P) THESE FILES ARE IN PSYCHRO4.NMF */ /* ENTHSAP (TDB,P) */ /* RHODRYP (TDB,W,P) */ /* RHOMOIP (TDB,W,P) */ /* */ /* DEWPNTPJ */ /* ENTHSAPJ */ /* RHODRYPJ */ /* RHOMOIPJ */ /* */ /*********************************************************************/ FUNCTION FLOAT DEWPNTP (W,P) /*********************************************************************/ /* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations */ /*********************************************************************/ /* FUNCTION: DEWPNT */ /* */ /* LANGUAGE: FORTRAN 77 */ /* */ /* PURPOSE: Calculate the dewpoint temperature given */ /* humidity ratio */ /* CODE IS TAKEN FROM ASHRAE HVAC 2 TOOLKIT */ /*********************************************************************/ /* INPUT VARIABLES: */ /* W Humidity ratio (-) */ /* Patm Atmospheric pressure (Pa) */ /* */ /* OUTPUT VARIABLES: */ /* DewPoint Dew point temperature of air (C) */ /* */ /* PROPERTIES: */ /*********************************************************************/ /* MAJOR RESTRICTIONS: None */ /* */ /* WRITTEN: Michael J. Brandemuehl */ /* */ /* ADDAPTED TO NMF BY: Mika Vuolle */ /* */ /* */ /* DATE: January 1, 1992 */ /* August 30, 1999 rev 1.0 */ /* */ /* INCLUDE FILES: globfor.inc */ /* SUBROUTINES CALLED: None */ /* FUNCTIONS CALLED: SATTEMP: Call SATPRES */ /* XITERATE */ /* REVISION HISTORY: */ /* */ /* */ /* REFERENCE: 1989 ASHRAE Handbook - Fundamentals */ /* */ /* ASHRAE HVAC 2 TOOLKIT */ /* Pages 7-3 - 7-4 */ /*********************************************************************/ /* INTERNAL VARIABLES: */ /* pw Partial water vapor pressure (Pa) */ /* small Small number */ /*********************************************************************/ LANGUAGE F77 INPUT FLOAT W,P; CODE DOUBLE PRECISION FUNCTION DewPntP (W,P) DOUBLE PRECISION W,pw,P,small DOUBLE PRECISION SATTEMP C INCLUDE 'globfor.inc' DATA small/1.E-9/ C1*** Test for "dry" air IF (W .LT. small) THEN DewPntP = -999 ELSE C1*** Calculate the partial water vapor pressure as a function of C1*** humidity ratio. pw= P*W/(.62198+W) C1*** Calculate dewpoint as saturation temperature at water vapor C1*** partial pressure DewPntP = SATTEMP(pw) ENDIF 999 RETURN END END_CODE EXTENSIONS JACOBIAN DEWPNTPJ; JACOBIAN_UNDER W,P; END_EXTENSIONS END_FUNCTION FUNCTION VOID DEWPNTPJ (W,P,J) LANGUAGE F77 INPUT FLOAT W,P; OUTPUT FLOAT J; CODE SUBROUTINE DEWPNTPJ(W,P,J) C INCLUDE 'globfor.inc' DOUBLE PRECISION W,P,J(2) DOUBLE PRECISION dW,DEWPNTP, small,dP DATA small/1.E-9/ IF (W == 0) THEN dW = .01 ELSE dW = 0.01*W ENDIF dP = 0.01*P IF (W .LE. SMALL) THEN J(1) = 0.0 J(2) = 0.5 * (DEWPNTP(SMALL,P+dP) - DEWPNTP(SMALL,P-dP))/dP ELSE J(1) = 0.5 * (DEWPNTP(W+dW,P) - DEWPNTP(W-dW,P))/dW J(2) = 0.5 * (DEWPNTP(W,P+dP) - DEWPNTP(W,P-dP))/dP ENDIF RETURN END END_CODE END_FUNCTION FUNCTION FLOAT ENTHSAP (TDB,P) /*********************************************************************/ /* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations */ /*********************************************************************/ /* FUNCTION: ENTHSAT */ /* */ /* LANGUAGE: FORTRAN 77 */ /* */ /* PURPOSE: Calculate the enthalpy at saturation for */ /* given dry bulb temperature. */ /* Atmospheric pressure is NOT constant */ /* CODE IS TAKEN FROM ASHRAE HVAC 2 TOOLKIT */ /*********************************************************************/ /* INPUT VARIABLES: */ /* TDB Dry bulb temperature (C) */ /* P Atmospheric pressure (Pa) */ /* */ /* OUTPUT VARIABLES: */ /* EnthSat Enthalpy at saturation (J/kg) */ /* */ /* PROPERTIES: */ /*********************************************************************/ /* MAJOR RESTRICTIONS: None */ /* */ /* WRITTEN: Shauna Gabel */ /* Michael J. Brandemuehl */ /* */ /* ADDAPTED TO NMF BY: Mika Vuolle */ /* */ /* */ /* DATE: January 1, 1992 */ /* August 31, 1999 rev 1.0 */ /* */ /* INCLUDE FILES: None */ /* SUBROUTINES CALLED: None */ /* FUNCTIONS CALLED: SATPRES */ /* HUMRAT */ /* ENTHAL */ /* */ /* REVISION HISTORY: None */ /* */ /* REFERENCE: 1989 ASHRAE Handbook - Fundamentals */ /* */ /* ASHRAE HVAC 2 TOOLKIT */ /* Pages 7-9 - 7-10 */ /*********************************************************************/ /* INTERNAL VARIABLES: */ /* psat Saturated water vapor pressure (Pa) */ /* w Humidity ratio (-) */ /*********************************************************************/ LANGUAGE F77 INPUT FLOAT TDB; CODE DOUBLE PRECISION FUNCTION ENTHSAP (TDB,P) DOUBLE PRECISION SATPRES,HUMRAT,ENTHAL DOUBLE PRECISION TDB,psat,w,P C1*** Calculate the saturation pressure at the given temperature. psat = SATPRES (TDB) C1*** Calculate the humidity ratio from the saturation pressure w = HUMRAT (P,psat) C1*** Calculate the enthalpy as a function of dry bulb temperature C1*** and humidity ratio. ENTHSAP = ENTHAL (TDB,w) RETURN END END_CODE EXTENSIONS JACOBIAN ENTHSAPJ; JACOBIAN_UNDER TDB,P; END_EXTENSIONS END_FUNCTION FUNCTION VOID ENTHSAPJ (TDB,P,J) LANGUAGE F77 INPUT FLOAT TDB,P; OUTPUT FLOAT J; CODE SUBROUTINE ENTHSAPJ (TDB,P,J) DOUBLE PRECISION TDB,P,J(2),EnthsaP DOUBLE PRECISION dTDB, dP IF (TDB == 0) THEN dTDB = .01 ELSE dTDB = 0.01*TDB ENDIF dP = 0.01*P J(1) = 0.5*(EnthsaP(TDB+dTDB,P)-EnthsaP(TDB-dTDB,P))/dTDB J(2) = 0.5*(EnthsaP(TDB,P+dP)-EnthsaP(TDB,P-dP))/dP RETURN END END_CODE END_FUNCTION FUNCTION FLOAT RHODRYP (TDB,W,P) /************************************************************************/ /* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations */ /************************************************************************/ /* FUNCTION: RHODRY */ /* */ /* LANGUAGE: FORTRAN 77 */ /* */ /* PURPOSE: Calculate dry air density. */ /* Atmospheric pressure is NOT constant */ /* */ /************************************************************************/ /* INPUT VARIABLES */ /* TDB Dry bulb temperature (C)*/ /* W Humidity ratio (-)*/ /* P Atmospheric pressure (Pa)*/ /* */ /* OUTPUT VARIABLES */ /* RhoDry Density of dry air (kg/m3)*/ /* */ /* PROPERTIES */ /* RAir Gas constant for air (J/kg C)*/ /************************************************************************/ /* MAJOR RESTRICTIONS: Perfect gas relationships */ /* */ /* DEVELOPER: Shauna Gabel */ /* Michael J. Brandemuehl, PhD, PE */ /* University of Colorado at Boulder */ /* */ /* ADDAPTED TO NMF BY: Mika Vuolle */ /* */ /* */ /* DATE: January 1, 1992 */ /* April 31, 1999 rev 1.0 */ /* */ /* */ /* INCLUDE FILES: globfor.inc */ /* SUBROUTINES CALLED: None */ /* FUNCTIONS CALLED: None */ /* */ /* REVISION HISTORY: */ /* */ /* REFERENCE: 1989 ASHRAE Handbook - Fundamentals */ /* */ /* ASHRAE HVAC 2 TOOLKIT */ /* Pages 7-17 - 7-18 */ /************************************************************************/ /* INTERNAL VARIABLES: */ /* pAir Partial pressure of dry air (Pa)*/ /************************************************************************/ LANGUAGE F77 INPUT FLOAT TDB,W,P; CODE DOUBLE PRECISION FUNCTION RHODRYP (TDB,W,P) DOUBLE PRECISION TDB,W,pAir,P INCLUDE 'globfor.inc' C1*** Calculate the dry air density from perfect gas laws. IF (W .LT. 0) THEN pAir = P ELSE pAir = 0.62198*P/(0.62198+W) ENDIF IF (TDB .LT. -200) THEN RhoDryP = p/RAir/(-200-ABS_ZERO) ELSE RhoDryP = p/RAir/(TDB-ABS_ZERO) ENDIF RETURN END END_CODE EXTENSIONS JACOBIAN RHODRYPJ; JACOBIAN_UNDER TDB,W,P; END_EXTENSIONS END_FUNCTION FUNCTION VOID RHODRYPJ (TDB,W,P,J) LANGUAGE F77 INPUT FLOAT TDB,W,P; OUTPUT FLOAT J; CODE SUBROUTINE RHODRYPJ (TDB,W,P,J) DOUBLE PRECISION TDB,W,P,J(3) INCLUDE 'globfor.inc' IF ((W .LT. 0).AND. (TDB .LT. -200)) THEN J(1) = 0 J(2) = 0 J(3) = 1.0/(RAir*(-200-ABS_ZERO)) ELSEIF (W .LT. 0) THEN J(1) = -Patm/(RAir*(TDB-ABS_ZERO)**2) J(2) = 0 J(3) = 1.0/(RAir*(TDB-ABS_ZERO)) ELSEIF (TDB .LT. -200) THEN J(1) = 0.0 J(2) = -0.62198*Patm/(RAir*(-200-ABS_ZERO)*(0.62198+W)**2) J(3) = 0.62198/(RAir*(-200-ABS_ZERO)*(0.62198+W)) ELSE J(1) = -0.62198*Patm/((0.62198+W)*RAir*(TDB-ABS_ZERO)**2) J(2) = -0.62198*Patm/(RAir*(TDB-ABS_ZERO)*(0.62198+W)**2) J(3) = 0.62198/(RAir*(TDB-ABS_ZERO)*(0.62198+W)) ENDIF RETURN END END_CODE END_FUNCTION FUNCTION FLOAT RHOMOIP (TDB,W,P) /************************************************************************/ /* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations */ /************************************************************************/ /* FUNCTION: RHOMOIO */ /* */ /* LANGUAGE: FORTRAN 77 */ /* */ /* PURPOSE: Calculate moist air density from */ /* dry air density and humidity ratio */ /* Atmospheric pressure is NOT constant */ /************************************************************************/ /* INPUT VARIABLES */ /* TDB Dry bulb temperature (C)*/ /* W Humidity ratio (-)*/ /* Patm Atmospheric pressure (Pa)*/ /* */ /* OUTPUT VARIABLES */ /* RhoMois Density of moist air (kg/m3)*/ /* */ /* PROPERTIES */ /* RAir Gas constant for air (J/kg C)*/ /************************************************************************/ /* MAJOR RESTRICTIONS: None */ /* */ /* DEVELOPER: Shauna Gabel */ /* Michael J. Brandemuehl, PhD, PE */ /* University of Colorado at Boulder */ /* */ /* ADDAPTED TO NMF BY: Mika Vuolle */ /* */ /* */ /* DATE: August 25, 1999 rev 1.0 */ /* */ /* */ /* SUBROUTINES CALLED: None */ /* FUNCTIONS CALLED: None */ /* */ /* REVISION HISTORY: None */ /* */ /* */ /* REFERENCE: 1989 ASHRAE Handbook - Fundamentals */ /* */ /* ASHRAE HVAC 2 TOOLKIT */ /* Pages 7-19 - 7-20 */ /************************************************************************/ /* INTERNAL VARIABLES: */ /* pAir Partial pressure of dry air (Pa)*/ /************************************************************************/ LANGUAGE F77 INPUT FLOAT TDB,W,P; CODE DOUBLE PRECISION FUNCTION RHOMOIP (TDB,W,P) DOUBLE PRECISION TDB,W,pAir,P INCLUDE 'globfor.inc' C1*** Calculate the moist air density. IF (W .LE. 0) THEN pAir = P ELSE pAir = 0.62198*P/(0.62198+W) ENDIF IF ((TDB .LT. -200) .AND. (W .LE. 0))THEN RhoMoiP = pAir/RAir/(-200-ABS_ZERO) ELSEIF (TDB .LT. -200) THEN RhoMoiP = pAir/RAir/(-200-ABS_ZERO)*(1.+W) ELSEIF (W .LE. 0.0) THEN RhoMoiP = pAir/RAir/(TDB-ABS_ZERO) ELSE RhoMoiP = pAir/RAir/(TDB-ABS_ZERO)*(1.+W) ENDIF RETURN END END_CODE EXTENSIONS JACOBIAN RHOMOIPJ; JACOBIAN_UNDER TDB,W,P; END_EXTENSIONS END_FUNCTION FUNCTION VOID RHOMOIPJ (TDB,W,P,J) LANGUAGE F77 INPUT FLOAT TDB,W,P; OUTPUT FLOAT J; CODE SUBROUTINE RHOMOIPJ (TDB,W,P,J) DOUBLE PRECISION TDB,W,P,J(3) INCLUDE 'globfor.inc' IF ((W .LE. 0) .AND. (TDB .LT. -200)) THEN J(1) = 0 J(2) = 0 J(3) = 1./(RAir*(-200-ABS_ZERO)) ELSEIF (TDB .LT. -200) THEN J(1) = 0 J(2) = 0.62198*P/(RAir*(-200-ABS_ZERO)*(0.62198+W)) & -0.62198*P*(1.+W)/ & (RAir*(-200-ABS_ZERO)*(0.62198+W)**2) J(3) = 0.62198*(1.+W)/(RAir*(-200-ABS_ZERO)*(0.62198+W)) ELSEIF (W .LE. 0) THEN J(1) = -P/(RAir*(TDB-ABS_ZERO)**2) J(2) = 0 J(3) = 1./(RAir*(TDB-ABS_ZERO)) ELSE J(1) = -0.62198*P*(1.+W)/((0.62198+W)*RAir*(TDB-ABS_ZERO)**2) J(2) = 0.62198*P/(RAir*(TDB-ABS_ZERO)*(0.62198+W)) & -0.62198*P*(1.+W)/(RAir*(TDB-ABS_ZERO)*(0.62198+W)**2) J(3) = 0.62198*(1.+w)/(RAir*(TDB-ABS_ZERO)*(0.62198+W)) ENDIF RETURN END END_CODE END_FUNCTION FUNCTION FLOAT WETBULP (TDB,W,P) /**************************************************************************/ /* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations */ /**************************************************************************/ /* FUNCTION: WETBULB */ /* */ /* LANGUAGE: FORTRAN 77 */ /* */ /* PURPOSE: Calculate wet bulb temperature from dry */ /* bulb temperature and humidity ratio */ /* Pressure is NOT constant */ /* */ /**************************************************************************/ /* INPUT VARIABLES */ /* TDB Dry bulb temperature (C) */ /* W Humidity ratio of air (-) */ /* P Atmospheric pressure (Pa) */ /* */ /* OUTPUT VARIABLES */ /* WetBulb Wet bulb temperature (C) */ /* */ /* PROPERTIES: */ /* Hfg Latent heat of vaporization of water (J/kg) */ /* CpAir Specific heat of air (J/kg C) */ /* CpVap Specific heat of water vapor (J/kg C) */ /* CpWat Specific heat of water (J/kg C) */ /**************************************************************************/ /* MAJOR RESTRICTIONS: None */ /* */ /* DEVELOPER: Shauna Gabel */ /* Michael J. Brandemuehl, PhD, PE */ /* University of Colorado at Boulder */ /* */ /* ADDAPTED TO NMF BY: Mika Vuolle */ /* */ /* */ /* DATE: January 1, 1992 */ /* August 30, 1999 */ /* */ /* INCLUDE FILES: globfor.inc */ /* SUBROUTINES CALLED: None */ /* FUNCTIONS CALLED: SATPRES */ /* HUMRATIO */ /* SATTEMP */ /* XITERATE */ /* */ /* REVISION HISTORY: None */ /* */ /* REFERENCE: 1989 ASHRAE Handbook - Fundamentals */ /**************************************************************************/ /* INTERNAL VARIABLES: */ /* tBoil Boiling temperature of water at given pressure (C) */ /* psatStar Saturation pressure at wet bulb temperature (C) */ /* wStar Humidity ratio as a function of PsatStar (-) */ /* newW Humidity ratio calculated with wet bulb guess (-) */ /* error Deviation of dependent variable in iteration */ /* iter Iteration counter */ /* icvg Iteration convergence flag */ /* F1,F2 Previous values of dependent variable in XITERATE */ /* X1,X2 Previous values of independent variable in XITERATE */ /**************************************************************************/ LANGUAGE F77 INPUT FLOAT TDB,W,P; CODE DOUBLE PRECISION FUNCTION WETBULP (TDB,W,P) DOUBLE PRECISION TDB,W,P,newW DOUBLE PRECISION tBoil,psatStar,wStar,error DOUBLE PRECISION SATTEMP, SATPRES, HUMRAT, XITERATE DOUBLE PRECISION X1,F1,X2,F2 INTEGER iter, itmax, icvg DATA itmax/20/ INCLUDE 'globfor.inc' C1*** Initial temperature guess tBoil = SATTEMP (P) WetBulP = MAX( MIN(WetBulP,TDB,(tBoil-0.1)), 0.) C1*** Begin iteration loop DO 100 iter = 1,itmax IF (WetBulP .GE. (tBoil-0.09) ) WETBULP = tBoil-0.1 C1*** Determine the saturation pressure for wet bulb temperature psatStar = SATPRES (WetBulP) C1*** Determine humidity ratio for given saturation pressure wStar = HUMRAT (P,psatStar) C1*** Calculate new humidity ratio and determine difference from known C1*** humidity ratio newW = ((Hfg-(CpWat-CpVap)*WetBulP)*wStar- & CpAir*(TDB-WetBulP))/(Hfg+CpVap*TDB & -CpWat*WetBulP) C1*** Check error, if not satisfied, calculate new guess and iterate error = W-newW WetBulP = XITERATE(WetBulP,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 C PRINT OUT REMOVED BY MIKA VUOLLE C 980616 C WRITE(*,1009) itmax C1009 FORMAT(/1X,'*** ERROR IN FUNCTION WetBulb ***'/ C & 1X,' Wet bulb temperature has not ' C & 'converged after ',I2,' iterations'/) 900 IF (WetBulP .GT. TDB) WetBulP = TDB 999 RETURN END END_CODE EXTENSIONS JACOBIAN WETBULPJ; JACOBIAN_UNDER TDB,W,P; END_EXTENSIONS END_FUNCTION FUNCTION VOID WETBULPJ(TDB,W,P,J) LANGUAGE F77 INPUT FLOAT TDB,W,P; OUTPUT FLOAT J; CODE SUBROUTINE WETBULPJ(TDB,W,P,J) DOUBLE PRECISION TDB,W,P,J(3),WetBulP DOUBLE PRECISION dTDB,dW,dP IF (TDB == 0) THEN dTDB=.01 ELSE dTDB=0.01*TDB ENDIF IF (W == 0) THEN dW=0.01 ELSE dW=0.01*W ENDIF dP=0.01*P J(1)=0.5*(WetBulP(TDB+dTDB,W,P)-WetBulP(TDB-dTDB,W,P))/dTDB J(2)=0.5*(WetBulP(TDB,W+dw,P)-WetBulP(TDB,W-dw,P))/dw J(3)=0.5*(WetBulP(TDB,W,P+dP)-WetBulP(TDB,W,P-dP))/dP RETURN END END_CODE END_FUNCTION