ALGORITHMIC_MODEL SyntClim ABSTRACT "Generate SYNThetic data for input into CLIMate model. Purpose: Generate synthetic data for climate model. The model generates the same data as a climate file; i.e. it calculates and delivers the following data: PAir Atmospheric pressure - 1E5 TAir Air temperature RelHum Relative humidity WindDir Direction of wind WindVel Velocity of wind IDirNorm Direct normal radiation IDiffHor Diffuse radiation on a horizontal surface Humidity can be specified in either of two ways 1. Minimum relative humidity (occurring at maximum dry bulb temperature) can be given. 2. For summer conditions, desired design TWetBulb can be specified. Values can be fetched e.g from ASHRAE Fundamentals or CIBS Guide. Absolute humidity is kept constant, except when the specified drybulb variation leads to saturation. In this case, the relative humidity is set to 100% and absolute humidity is assumed to be restored as fast as rising temperature allows. For documentation on integration time, cf CLIMATE.NMF. Limits: Direct radiation is presented on horizontal surface. Date: January 24, 1998 Made by: Axel Bring Call: EnthSat, HumRat, HumTH, SatPres Reference: 'CIBS Guide, Design Temperatures - Approximate Method, pp A2.18-20', 'ASHRAE Fundamentals, Weather Data, chapter 24'. Changes: 980901 MV A and B1 local variables for ASHRAE method calc are added. 980921 MV Equations for BRIS winter case are added. Diffuse sky factor of ASHRAE method is added as a function of month Reduction factor is added for direct normal radiation 991111 AB (MV) The wetbulb calc is updated. Old version generates too high moisture contents. It used wetbulb as a dewpoint. " /* Revision history AB 001031 Avoid overflow in EXP when ElevSun < 0.1 Deg. AB 980124 1st draft */ EQUATIONS /* ***************** Copy from CLIMATE.NMF ******************** */ /* CALCULATE SUN POSITION */ /* Transform from time in Universal Time to hours within year */ /* Calculate day number and hour, starting from January 1st */ /* Remove four year periods, after the non leap year 1900; then remove full 365 day years. The last day in a leap year will wrap around to the 1st. */ Seconds := Time; DayFrom1901 := Seconds/24/3600 - 365; If DayFrom1901 < 0 THEN /* year 0 */ DayInYear := DayFrom1901 + 365; ELSE /* year 1.., calulate leap year; can be Universal Time for ICE */ DayInFourYears := DayFrom1901 - Int(DayFrom1901/1461) * 1461; DayInYear := DayInFourYears - Int(DayInFourYears/365) * 365; END_IF; IF DayInYear < 0 OR DayInYear > 365 THEN CALL nmf_error ("Integration time must be non-negative"); END_IF; TimeHr := DayInYear * 24; DayNr := INT (DayInYear) + 1; Hour := TimeHr - 24*(DayNr-1); /* Help variables to calculate solar time */ Declination := 23.45*deg2rad * SIN(2*PI * (284+DayNr)/365); B := 2*PI * (DayNr-81)/364; E := 9.87*SIN(2*B) - 7.53*COS(B) - 1.5*SIN(B); SolTime := Hour + 4*(LongTimezone - Long)/60 + E/60; /* Check if SolTime is <0 or >24 */ SolarTime := IF (SolTime < 0) THEN SolTime + 24; ELSE_IF (SolTime > 24) THEN SolTime - 24; ELSE SolTime END_IF; /* Calculate hour angle from south, positive westward */ Omega := (SolarTime-12)*15*deg2rad; /* Elevation angle of Sun */ /* ASHRAE Fundamentals Chapter 27 */ SinElevSun := COS(Lat*deg2rad)*COS(Declination)*COS(Omega) + SIN(Lat*deg2rad)*SIN(Declination); ElevSun := ASIN(SinElevSun)*rad2deg; /* ************** End of copy from CLIMATE.NMF ***************** */ PAir := pAirGiven; XAir := xAirGiven; /* Calculate temperatures and humidity ratio */ TAir := tDryMean + tDryAmpl * sin ((SolarTime - hrMaxTemp)*15*deg2rad + PI/2); /* Saturation pressure at TAir */ PSat := SatPres (TAir); RelHumTry := IF relHumAtMax < 0.01 THEN /* Use wetbulb temperature */ humRatio / HumRat (1E5+PAir, Psat) * 100 ELSE /* Use relative humidity at max drybulb */ PVapMax / PSat * 100 END_IF; RelHum := IF RelHumTry < 100 THEN RelHumTry ELSE 100 END_IF; /* Average wind conditions */ WindDir := dominDir; WindVelRef := avgWindVel; /* Factors of ASHRAE radiation calculation */ IF DayNr < 31 THEN A := 1230; B1 := 0.142; C := 0.058; ELSE_IF DayNr < 59 THEN A := 1215; B1 := 0.144; C := 0.060; ELSE_IF DayNr < 90 THEN A := 1186; B1 := 0.156; C := 0.071; ELSE_IF DayNr < 120 THEN A := 1136; B1 := 0.180; C := 0.097; ELSE_IF DayNr < 151 THEN A := 1104; B1 := 0.196; C := 0.121; ELSE_IF DayNr < 181 THEN A := 1088; B1 := 0.205; C := 0.134; ELSE_IF DayNr < 212 THEN A := 1085; B1 := 0.207; C := 0.136; ELSE_IF DayNr < 243 THEN A := 1107; B1 := 0.201; C := 0.122; ELSE_IF DayNr < 273 THEN A := 1151; B1 := 0.177; C := 0.092; ELSE_IF DayNr < 304 THEN A := 1192; B1 := 0.160; C := 0.073; ELSE_IF DayNr < 334 THEN A := 1221; B1 := 0.149; C := 0.063; ELSE A := 1223; B1 := 0.142; C := 0.057; END_IF; /* Calculate solar radiation */ IDirNorm := RedFac * IF rad_model < 0.5 THEN /* Solar intensity by ASHRAE Fundamentals Chapter 27 */ /* avoid overflow in EXP when ElevSun < 0.1 Deg */ IF ElevSun < 0.1 then 0 ELSE A/EXP(B1/SIN(ElevSun*Pi/180)) END_IF; ELSE IF DayNr < 120 OR DayNr > 273 THEN /* Solar intensity in winter by BRIS */ IF ElevSun < 0 then 0 ELSE IF ElevSun < 15 THEN 1.163*ElevSun*(100.766+ElevSun*(-8.1631+ElevSun* (0.38131+ElevSun*(-0.007066)))) ELSE 1.163*921/Exp(.109/sin(ElevSun*Pi/180)) END_IF END_IF /* Solar intensity in summer by BRIS */ ELSE IF ElevSun < 0 then 0 ELSE IF ElevSun < 15 THEN 1.163*ElevSun*(87.616+ElevSun*(-6.9947+ElevSun* (0.32331+ElevSun*(-0.005799)))) ELSE 1.163*921/Exp(.139/sin(ElevSun*Pi/180)) END_IF END_IF END_IF; END_IF; IDiffuse := IF rad_model < 0.5 THEN /* ASHRAE model */ C*IDirNorm ELSE IF ElevSun < 60 THEN -0.82189+ElevSun*(5.26273+ ElevSun*(-0.09372+ElevSun*0.0006)) ELSE (ElevSun-60) / 30 * (110-107.15)+107.15 END_IF; END_IF; IDiffHor := IF IDiffuse < 0 THEN 0 ELSE IDiffuse END_IF; LINKS /* type name variables .... */ /* Output data, equivalent to climate file */ ClimData DataOut PAir, TAir, XAir, RelHum, WindDir, WindVelRef, IDirNorm, IDiffHor; VARIABLES /* type name role def min max description*/ Pressure PAir OUT 1325 SMALL BIG "Atmospheric pressure" Temp TAir OUT 20 ABS_ZERO BIG "Temperature of air" Fraction_y XAir OUT 594 0 BIG "CO2 fraction" Factor RelHum OUT 50 0 100 "Rel humidity of air" Factor RelHumTry LOC 50 0 100 "Tentative rel humidity of air" Angle WindDir OUT 0 0 360 "Direction of wind" Vel WindVelRef OUT 1 0 BIG "Speed of meteorological wind" RadA IDirNorm OUT 0 0 BIG "Direct normal rad" RadA IDiffuse LOC 0 0 BIG "1st diffuse rad" RadA IDiffHor OUT 0 0 BIG "Diffuse rad on hor surf" Angle ElevSun LOC 27.0 -90 90 "Elevation angle of sun" Pressure PSat LOC 2365 SMALL BIG "Saturation pressure of water vapor" Generic DayNr LOC 80 1 366 "Whole day number in year, 1 = January 1st" Generic DayInYear LOC 80 1 366 "Broken day nr in year, 0 = January 1st" Generic DayInFourYears LOC 80 1 366 "Broken day nr in four year leap period, 0 = January 1st" Generic DayFrom1901 LOC 80 1 366 "Broken day nr from 1901-01-01, 0-.." Generic Hour LOC 14 0 24 "Hour number in day, std clock time" Generic Seconds LOC 86400 0 BIG "Time from 1900-01-01" Generic TimeHr LOC 0 -BIG BIG "Integration time [hr]" AngleR Declination LOC -0.007 SMALL BIG "Declination angle" Angle B LOC -0.017 -BIG BIG "Help var for solar time" Angle E LOC -7.84 -BIG BIG "Help var for solar time" AngleR Omega LOC 0.402 -BIG BIG "Hour angle from south, 15 deg/hour" Hour SolTime LOC 13.54 0 24 "Help var for solar time" Hour SolarTime LOC 13.54 0 24 "Local solar time" Factor SinElevSun LOC 0.454 -1 1 "Help var, SIN(ElevSun)" RadA A LOC 1088 1085 1233 "Factor of ASHRAE method" Factor B1 LOC 0.205 0.142 0.207 "Factor of ASHRAE method" Factor C LOC 0.136 0.057 0.136 "Sky diffuse factor of ASHRAE method " PARAMETERS /* type name role def min max description*/ Pressure pAirGiven S_P 1325 -50000 BIG "(Atmospheric press)-1E5" Temp tDryMean S_P 20 ABS_ZERO BIG "Mean temperature of air" Temp tDryAmpl S_P 5 0 50 "Amplitude of air temp " Generic hrMaxTemp S_P 15 0 24 "Hour (solar) when temp peaks" Fraction_y xAirGiven S_P 594 0 BIG "CO2 fraction" Factor relHumAtMax S_P 50 0 100 "Rel humiditiy at max TDry Specify 0 if tWetBulb given" Temp tWetBulb S_P 20 ABS_ZERO BIG "Design wet bulb temp Only used if relHumAtMax = 0" Angle dominDir S_P 0 0 360 "Dominating direction of wind" Vel avgWindVel S_P 2 0 BIG "Avg speed of meteorological wind" Angle Lat S_P 60 -90 90 "Local latitude" Angle Long S_P -25 -180 180 "Local longitude" Angle LongTimezone S_P -30 -180 180 "Time zone longitude" Factor diff_model S_P 0 0 1 "0 = ASHRAE 1 = Threlkeld" Factor rad_model S_P 0 0 1 "0 = ASHRAE 1 = BRIS" Factor RedFac S_P 1 0 1.25 "Clearness number (ie Reduction factor)" Factor deg2rad C_P 0.0175 SMALL BIG "Conversion factor Deg to Rad" Factor rad2deg C_P 57.296 SMALL BIG "Conversion factor Rad to Deg" /* Calculated parameters */ Temp tDryMax C_P "Max drybulb temp" Factor humRatio C_P "Humidity ratio, used except at saturation" Pressure pSatMax C_P 1182 SMALL BIG "Saturation pressure of water vapor at max temp" Pressure pVapMax C_P 1182 SMALL BIG "Partial pressure of water vapor at max temp" PARAMETER_PROCESSING /* Conversion factor from Deg to Rad and vice-versa */ deg2rad := PI/180.0; rad2deg := 180.0/PI; tDryMax := tDryMean + tDryAmpl; /* Vapour pressure */ pSatMax := SatPres (tDryMax); /* Select mode to specify humidity */ IF relHumAtMax < 0.01 THEN /* Use wetbulb temperature */ humRatio := HumTH (tDryMax, EnthSat (tWetBulb)); ELSE /* Use relative humidity at max drybulb */ pVapMax := relHumAtMax/100 * PsatMax; END_IF; END_MODEL