WSL/SLF GitLab Repository

Atmosphere.cc 28.2 KB
 Mathias Bavay committed Jan 26, 2011 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 ``````/***********************************************************************************/ /* Copyright 2009 WSL Institute for Snow and Avalanche Research SLF-DAVOS */ /***********************************************************************************/ /* This file is part of MeteoIO. MeteoIO is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. MeteoIO is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with MeteoIO. If not, see . */ #include #include #include #include `````` 23 ``````#include `````` Mathias Bavay committed Jan 26, 2011 24 25 26 ``````#include namespace mio { `````` 27 28 29 30 ``````/** * @brief Calculate the black body emissivity * @param lwr longwave radiation emitted by the body (W m-2) * @param T surface temperature of the body (K) `````` Mathias Bavay committed Nov 07, 2011 31 `````` * @return black body emissivity (0-1) `````` 32 33 34 `````` */ double Atmosphere::blkBody_Emissivity(const double& lwr, const double& T) { const double T2 = T*T; `````` 35 36 37 38 39 `````` const double ea = lwr / (Cst::stefan_boltzmann * (T2*T2)); if(ea > 1.0) return 1.0; return ea; `````` 40 41 42 43 44 45 ``````} /** * @brief Calculates the black body long wave radiation knowing its emissivity * @param ea emissivity of the body (0-1) * @param T surface temperature of the body (K) `````` Mathias Bavay committed Nov 07, 2011 46 `````` * @return black body radiation (W/m^2) `````` 47 48 49 50 51 52 `````` */ double Atmosphere::blkBody_Radiation(const double& ea, const double& T) { const double T2 = T*T; return ( ea * (Cst::stefan_boltzmann * (T2*T2)) ); } `````` Mathias Bavay committed Jan 26, 2011 53 54 55 56 57 58 59 60 61 ``````/** * @brief Standard atmosphere pressure * @param altitude altitude above sea level (m) * @return standard pressure (Pa) */ double Atmosphere::stdAirPressure(const double& altitude) { const double p0 = Cst::std_press; // Air and standard pressure in Pa const double lapse_rate = 0.0065; // K m-1 const double sea_level_temp = 288.15; // K `````` 62 `````` const double expo = Cst::gravity / (lapse_rate * Cst::gaz_constant_dry_air); `````` Mathias Bavay committed Jan 26, 2011 63 `````` const double R0 = Cst::earth_R0; // Earth's radius in m `````` Mathias Bavay committed Aug 16, 2011 64 `````` `````` Mathias Bavay committed Jan 26, 2011 65 `````` const double p = p0 * pow( 1. - ( (lapse_rate * R0 * altitude) / (sea_level_temp * (R0 + altitude)) ), expo ); `````` Mathias Bavay committed Aug 16, 2011 66 `````` `````` Mathias Bavay committed Jan 26, 2011 67 68 69 `````` return(p); } `````` 70 71 72 73 74 75 76 77 78 79 80 81 ``````/** * @brief Standard dry air pressure * @param altitude altitude above sea level (m) * @param temperature air temperature (K) * @return standard pressure (Pa) */ double Atmosphere::stdDryAirDensity(const double& altitude, const double& temperature) { return stdAirPressure(altitude)/(Cst::gaz_constant_dry_air*temperature); } /** * @brief Calculates the water vapor density, for a given temperature and vapor pressure `````` Mathias Bavay committed Mar 15, 2011 82 ``````* @param Temperature air temperature (K) `````` 83 84 85 86 87 88 89 ``````* @param VaporPressure water vapor pressure (Pa) * @return water vapor density (kg/m^3) */ double Atmosphere::waterVaporDensity(const double& Temperature, const double& VaporPressure) { return (Cst::water_molecular_mass*VaporPressure)/(Cst::gaz_constant*Temperature); } `````` 90 ``````/** `````` Mathias Bavay committed Aug 16, 2011 91 ``````* @brief Standard atmosphere wet bulb temperature. `````` 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 ``````* This gives the lowest temperature that could be reached by water evaporation. It is therefore linked to * relative humidity. This implementation assumes a standard atmosphere for pressure and saturation pressure. * @param T air temperature (K) * @param RH relative humidity (between 0 and 1) * @param altitude altitude above sea level (m) * @return wet bulb temperature (K) */ double Atmosphere::wetBulbTemperature(const double& T, const double& RH, const double& altitude) { const double L = Cst::l_water_vaporization; //latent heat of vaporisation const double mixing_ratio = Cst::gaz_constant_dry_air / Cst::gaz_constant_water_vapor; const double p = stdAirPressure(altitude); const double Vp = waterSaturationPressure(T); return ( T - (RH*Vp - Vp) * mixing_ratio * L / p / Cst::specific_heat_air ); } `````` Mathias Bavay committed Feb 14, 2012 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 ``````/** * @brief Wind log profile. * This is used to compute the equivalent wind speed at a different height. * @details It depends on the roughness length * that can take a default value (0.03) or any other, for example: *
* * * * * * * * * * * @param v_ref reference wind speed (m/s) * @param z_ref height of the reference wind speed (m) * @param z new height (m) * @param z0 roughness length (m) * @return wind speed at height z (m/s) */ double Atmosphere::windLogProfile(const double& v_ref, const double& z_ref, const double& z, const double& z0) { return v_ref * ( log(z/z0) / log(z_ref/z0) ); } `````` Mathias Bavay committed Feb 07, 2012 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 ``````/** * @brief Wind chill temperature. * This is an index aiming at expressing the human-percived feeling of air temperature on exposed skin due to wind. * This is NOT a sceintific measurement, only an index to express a subjective feeling. * It is inapplicable above 10 Celsius and a few m/s wind (5 m/s used here), therefore returning air temperature. * @param TA air temperature (K) * @param VW wind velocity (m/s) * @return wind chill temperature (K) */ double Atmosphere::windChill(const double& TA, const double& VW) { if(TA>(273.15+10.) || VW<5.) return TA; //not applicable in this case const double t = K_TO_C(TA); //in Celsius const double v = VW*3.6; //in km/h const double v016 = pow(v, 0.16); const double WCT = 13.12 + 0.6215*t - 11.37*v016 + 0.3965*t*v016; return C_TO_K(WCT); } /** * @brief Heat index. * This is an index aiming at expressing the human-perceived air temperature due to humidity. * This is NOT a sceintific measurement, only an index to express a subjective feeling. * It is inapplicable below 26.7 Celsius and below 40% humidity, therefore returning air temperature. * @param TA air temperature (K) * @param RH relative humidity (between 0 and 1) * @return Heat index (K) */ double Atmosphere::heatIndex(const double& TA, const double& RH) { if(TA<(273.15+26.7) || RH<0.4) return TA; //not applicable in this case const double t = K_TO_C(TA); //in Celsius const double t2 = t*t; const double rh = RH*100.; //in percent const double rh2 = rh*rh; const double c1=-8.784695, c2=1.61139411, c3=2.338549 , c4=-0.14611605; const double c5=-1.2308094e-2 , c6=-1.6424828e-2 , c7=2.211732e-3 , c8=7.2546e-4, c9=-3.582e-6; const double HI = c1 + c2*t + c3*rh + c4*t*rh + c5*t2 + c6*rh2 + c7*t2*rh + c8*t*rh2 + c9*t2*rh2; return C_TO_K(HI); } `````` Mathias Bavay committed Jan 26, 2011 184 185 186 187 188 189 190 191 ``````/** * @brief Standard water vapor saturation * @param T air temperature (K) * @return standard water vapor saturation pressure (Pa) */ double Atmosphere::waterSaturationPressure(const double& T) { double c2, c3; // varying constants `````` Mathias Bavay committed Jun 19, 2012 192 `````` if ( T < Cst::t_water_triple_pt ) { // for a flat ice surface `````` Mathias Bavay committed Jan 26, 2011 193 194 195 196 197 198 `````` c2 = 21.88; c3 = 7.66; } else { // for a flat water surface c2 = 17.27; c3 = 35.86; } `````` Mathias Bavay committed Jun 19, 2012 199 `````` const double exp_p_sat = c2 * (T - Cst::t_water_triple_pt) / (T - c3); //exponent `````` Mathias Bavay committed Jan 26, 2011 200 201 202 203 `````` return( Cst::p_water_triple_pt * exp( exp_p_sat ) ); } `````` Mathias Bavay committed Nov 07, 2011 204 205 206 207 208 209 210 211 212 213 214 215 216 ``````/** * @brief Virtual temperature multiplying factor. * In order to get a virtual temperature, multiply the original temperature by this factor. Note: * since e/p is actually used, both pressures only need to have the same units. * @param e vapor pressure (Pa) * @param p air pressure (Pa) * @return virtual temperature multiplying coefficient */ double Atmosphere::virtualTemperatureFactor(const double& e, const double& p) { const double epsilon = 0.622; return 1. / (1.-(1.-epsilon)*e/p); } `````` 217 ``````/** `````` 218 ``````* @brief Evaluate the atmosphere emissivity from the water vapor pressure and cloudiness. `````` 219 ``````* This is according to A. Omstedt, "A coupled one-dimensional sea ice-ocean model applied to a semi-enclosed basin", `````` Nander Wever committed Sep 24, 2013 220 ``````* Tellus, 42 A, 568-582, 1990, DOI:10.1034/j.1600-0870.1990.t01-3-00007.x `````` 221 222 ``````* @param RH relative humidity (between 0 and 1) * @param TA air temperature (K) `````` 223 ``````* @param cloudiness cloudiness (between 0 and 1, 0 being clear sky) `````` 224 ``````* @return emissivity (between 0 and 1) `````` 225 ``````*/ `````` 226 227 ``````double Atmosphere::Omstedt_emissivity(const double& RH, const double& TA, const double& cloudiness) { const double e0 = RH * waterSaturationPressure(TA); //water vapor pressure `````` 228 229 230 231 232 `````` const double eps_w = 0.97; const double a1 = 0.68; const double a2 = 0.0036; const double a3 = 0.18; `````` 233 `````` const double ea = (eps_w * (a1 + a2 * sqrt(e0)) * (1. + a3 * cloudiness * cloudiness)); //emissivity `````` Mathias Bavay committed Aug 16, 2011 234 `````` if(ea > 1.0) `````` 235 `````` return 1.0; `````` 236 237 238 239 `````` return ea; } /** `````` Mathias Bavay committed Aug 16, 2011 240 ``````* @brief Evaluate the long wave radiation from RH, TA and cloudiness. `````` 241 242 243 244 245 246 247 248 ``````* This is according to A. Omstedt, "A coupled one-dimensional sea ice-ocean model applied to a semi-enclosed basin", * Tellus, 42 A, 568-582, 1990, DOI:10.1034/j.1600-0870.1990.t01-3-00007. * @param RH relative humidity (between 0 and 1) * @param TA air temperature (K) * @param cloudiness cloudiness (between 0 and 1, 0 being clear sky) * @return long wave radiation (W/m^2) */ double Atmosphere::Omstedt_ilwr(const double& RH, const double& TA, const double& cloudiness) { `````` 249 `````` const double ea = Omstedt_emissivity(RH, TA, cloudiness); `````` 250 251 252 253 `````` return blkBody_Radiation(ea, TA); } /** `````` Mathias Bavay committed Aug 16, 2011 254 `````` * @brief Evaluate the atmosphere emissivity for clear sky. `````` 255 256 `````` * This uses the formula from Brutsaert -- "On a Derivable * Formula for Long-Wave Radiation From Clear Skies", Journal of Water Resources `````` 257 `````` * Research, 11, No. 5, October 1975, pp 742-744. `````` 258 `````` * Alternative: Satterlund (1979): Water Resources Research, 15, 1649-1650. `````` 259 `````` * @param RH relative humidity (between 0 and 1) `````` 260 261 262 `````` * @param TA Air temperature (K) * @return clear sky emissivity */ `````` 263 264 ``````double Atmosphere::Brutsaert_emissivity(const double& RH, const double& TA) { const double e0 = RH * waterSaturationPressure(TA); //water vapor pressure `````` 265 266 `````` const double e0_mBar = 0.01 * e0; const double exponent = 1./7.; `````` 267 `````` const double ea = 1.24 * pow( (e0_mBar / TA), exponent); `````` 268 `````` `````` 269 270 271 `````` if(ea>1.0) return 1.; return ea; `````` 272 273 274 ``````} /** `````` Mathias Bavay committed Aug 16, 2011 275 `````` * @brief Evaluate the long wave radiation for clear sky. `````` 276 277 `````` * This uses the formula from Brutsaert -- "On a Derivable * Formula for Long-Wave Radiation From Clear Skies", Journal of Water Resources `````` 278 `````` * Research, 11, No. 5, October 1975, pp 742-744. `````` 279 280 281 282 283 284 `````` * Alternative: Satterlund (1979): Water Resources Research, 15, 1649-1650. * @param RH relative humidity (between 0 and 1) * @param TA Air temperature (K) * @return long wave radiation (W/m^2) */ double Atmosphere::Brutsaert_ilwr(const double& RH, const double& TA) { `````` 285 `````` const double ea = Brutsaert_emissivity(RH, TA); `````` 286 `````` return blkBody_Radiation(ea, TA); `````` 287 288 ``````} `````` Mathias Bavay committed Mar 01, 2012 289 290 291 292 ``````/** * @brief Evaluate the atmosphere emissivity for clear sky. * This uses the formula from Dilley and O'Brien -- "Estimating downward clear sky * long-wave irradiance at the surface from screen temperature and precipitable water", `````` 293 `````` * Q. J. R. Meteorolo. Soc., 124, 1998, pp 1391-1401. The long wave is computed `````` Mathias Bavay committed Mar 01, 2012 294 295 296 297 298 299 300 301 `````` * and the ratio of this long wave to a black body emission gives an emissivity. * @param RH relative humidity (between 0 and 1) * @param TA near surface air temperature (K) * @return long wave radiation (W/m^2) */ double Atmosphere::Dilley_emissivity(const double& RH, const double& TA) { const double ilwr_dilley = Dilley_ilwr(RH, TA); const double ilwr_blkbody = blkBody_Radiation(1., TA); `````` 302 303 304 305 306 `````` const double ea = ilwr_dilley/ilwr_blkbody; if(ea>1.0) return 1.; return ea; `````` Mathias Bavay committed Mar 01, 2012 307 308 309 310 311 312 ``````} /** * @brief Evaluate the long wave radiation for clear sky. * This uses the formula from Dilley and O'Brien -- "Estimating downward clear sky * long-wave irradiance at the surface from screen temperature and precipitable water", `````` 313 `````` * Q. J. R. Meteorolo. Soc., 124, 1998, pp 1391-1401. `````` Mathias Bavay committed Mar 01, 2012 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 `````` * @param RH relative humidity (between 0 and 1) * @param TA near surface air temperature (K) * @return long wave radiation (W/m^2) */ double Atmosphere::Dilley_ilwr(const double& RH, const double& TA) { const double e0 = RH * waterSaturationPressure(TA) * 0.001; //water vapor pressure, kPa const double w = 4650.*e0/TA; //precipitable water, Prata 1996 const double tmp = TA/Cst::t_water_triple_pt; const double pow_tmp6 = tmp*tmp*tmp*tmp*tmp*tmp; return 59.38 + 113.7*pow_tmp6 + 96.96*sqrt(w/25.); } /** * @brief Evaluate the atmosphere emissivity for clear sky. * This uses the formula from Prata -- "A new long-wave formula for estimating `````` 330 `````` * downward clear-sky radiation at the surface", Q. J. R. Meteorolo. Soc., 122, 1996, pp 1127-1151. `````` Mathias Bavay committed Mar 01, 2012 331 332 333 334 335 336 337 `````` * @param RH relative humidity (between 0 and 1) * @param TA near surface air temperature (K) * @return long wave radiation (W/m^2) */ double Atmosphere::Prata_emissivity(const double& RH, const double& TA) { const double e0 = RH * waterSaturationPressure(TA) * 0.001; //water vapor pressure, kPa const double w = 4650.*e0/TA; //precipitable water, Prata 1996 `````` 338 `````` const double ea = 1. - (1.+w)*exp( -sqrt(1.2+3.*w) ); `````` Mathias Bavay committed Mar 01, 2012 339 `````` `````` 340 341 342 `````` if(ea>1.0) return 1.; return ea; `````` Mathias Bavay committed Mar 01, 2012 343 344 345 346 347 ``````} /** * @brief Evaluate the long wave radiation for clear sky. * This uses the formula from Prata -- "A new long-wave formula for estimating `````` 348 `````` * downward clear-sky radiation at the surface", Q. J. R. Meteorolo. Soc., 122, 1996, pp 1127-1151. `````` Mathias Bavay committed Mar 01, 2012 349 350 351 352 353 354 355 356 357 `````` * @param RH relative humidity (between 0 and 1) * @param TA near surface air temperature (K) * @return long wave radiation (W/m^2) */ double Atmosphere::Prata_ilwr(const double& RH, const double& TA) { const double epsilon = Prata_emissivity(RH, TA); return blkBody_Radiation(epsilon, TA); } `````` 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 ``````/** * @brief Evaluate the cloudiness from a given solar index. * This uses the formula from Kasten and Czeplak -- "Solar and terrestrial radiation * dependent on the amount and type of cloud", Sol. Energy, 24, 1980, pp 177-189. * The solar index is defined as measured radiation / clear sky radiation, values * outside of [0;1] will be truncated to [0;1]. * @param solarIndex solar index * @return cloudiness (between 0 and 1) */ double Atmosphere::Kasten_cloudiness(const double& solarIndex) { const double b1 = 0.75, b2 = 3.4; if(solarIndex>1.) return 0.; const double cloudiness = pow((1.-solarIndex)/b1, 1./b2); if(cloudiness>1.) return 1.; return cloudiness; } `````` 376 377 378 379 380 ``````/** * @brief Evaluate the long wave radiation for clear or cloudy sky. * This uses the formula from Crawford and Duchon -- "An Improved Parametrization * for Estimating Effective Atmospheric Emissivity for Use in Calculating Daytime * Downwelling Longwave Radiation", Journal of Applied Meteorology, `````` Mathias Bavay committed Oct 02, 2012 381 382 383 `````` * 38, 1999, pp 474-480. If no cloud cover fraction is provided, a parametrization * using iswr_meas and iswr_clear_sky will be used. These parameters can therefore safely * be set to IOUtils::nodata if cloudiness is provided. `````` 384 385 386 387 388 `````` * @param RH relative humidity (between 0 and 1) * @param TA Air temperature (K) * @param iswr_meas Measured Incoming Short Wave Radiation (W/m^2) * @param iswr_clear_sky Clear Sky Modelled Incoming Short Wave Radiation (W/m^2) * @param month current month (1-12, for a sinusoidal variation of the leading coefficients) `````` Mathias Bavay committed Oct 02, 2012 389 `````` * @param cloudiness Cloud cover fraction (between 0 and 1, optional) `````` 390 391 `````` * @return long wave radiation (W/m^2) or IOUtils::nodata at night time */ `````` Mathias Bavay committed Oct 02, 2012 392 393 394 395 396 397 398 399 400 401 402 403 404 ``````double Atmosphere::Crawford_ilwr(const double& RH, const double& TA, const double& iswr_meas, const double& iswr_clear_sky, const unsigned char& month, const double& cloudiness) { double clf; if(cloudiness==IOUtils::nodata) { if(iswr_meas<=0. || iswr_clear_sky<=0.) return IOUtils::nodata; clf = 1. - iswr_meas/iswr_clear_sky; //cloud fraction estimate if(clf<0.) clf=0.; } else { if(cloudiness<0. || cloudiness>1.) return IOUtils::nodata; clf = cloudiness; } `````` 405 406 407 408 409 `````` const double e = RH * waterSaturationPressure(TA); //near surface water vapor pressure const double e_mBar = 0.01 * e; const double epsilon = clf + (1.-clf) * (1.22 + 0.06*sin((month+2.)*Cst::PI/6.) ) * pow( (e_mBar/TA), 1./7.); `````` Mathias Bavay committed Mar 01, 2012 410 `````` return blkBody_Radiation(epsilon, TA); `````` 411 412 413 414 415 416 417 ``````} /** * @brief Evaluate the long wave radiation for clear or cloudy sky. * This uses the formula from Crawford and Duchon -- "An Improved Parametrization * for Estimating Effective Atmospheric Emissivity for Use in Calculating Daytime * Downwelling Longwave Radiation", Journal of Applied Meteorology, `````` Mathias Bavay committed Oct 02, 2012 418 419 420 `````` * 38, 1999, pp 474-480. If no cloud cover fraction is provided, a parametrization * using the current location (lat, lon, altitude) and ISWR will be used. These parameters can therefore safely * be set to IOUtils::nodata if cloudiness is provided. `````` 421 422 423 424 425 426 427 428 `````` * @param lat latitude of the point of observation * @param lon longitude of the point of observation * @param altitude altitude of the point of observation * @param julian julian date at the point of observation * @param TZ time zone at the point of observation * @param RH relative humidity (between 0 and 1) * @param TA Air temperature (K) * @param ISWR Measured Incoming Short Wave Radiation (W/m^2) `````` Mathias Bavay committed Oct 02, 2012 429 `````` * @param cloudiness Cloud cover fraction (between 0 and 1, optional) `````` 430 431 432 433 434 435 436 `````` * @return long wave radiation (W/m^2) or IOUtils::nodata at night time * Please note that this call might NOT be efficient for processing large amounts of points, * since it internally builds complex objects at every call. You might want to copy/paste * its code in order to process data in bulk. */ double Atmosphere::Crawford_ilwr(const double& lat, const double& lon, const double& altitude, const double& julian, const double& TZ, `````` Mathias Bavay committed Oct 02, 2012 437 `````` const double& RH, const double& TA, const double& ISWR, const double& cloudiness) `````` 438 ``````{ `````` Mathias Bavay committed Oct 02, 2012 439 `````` if(TA==IOUtils::nodata || RH==IOUtils::nodata) { `````` 440 441 442 443 444 445 446 `````` return IOUtils::nodata; } Date date(julian, TZ, 0.); int year, month, day; date.getDate(year, month, day); `````` Mathias Bavay committed Oct 02, 2012 447 448 449 450 451 `````` if(cloudiness==IOUtils::nodata) { if(ISWR==IOUtils::nodata) return IOUtils::nodata; SunObject Sun(lat, lon, altitude, julian, TZ); Sun.calculateRadiation(TA, RH, 0.5); //we force a terrain albedo of 0.5... `````` 452 453 454 `````` double toa, direct, diffuse; Sun.getHorizontalRadiation(toa, direct, diffuse); return Atmosphere::Crawford_ilwr(RH, TA, ISWR, direct+diffuse, static_cast(month)); `````` Mathias Bavay committed Oct 02, 2012 455 456 457 `````` } else { return Atmosphere::Crawford_ilwr(RH, TA, IOUtils::nodata, IOUtils::nodata, static_cast(month), cloudiness); } `````` 458 459 ``````} `````` 460 461 462 463 ``````/** * @brief Evaluate the long wave radiation for clear or cloudy sky. * This uses the formula from Unsworth and Monteith -- "Long-wave radiation at the ground", * Q. J. R. Meteorolo. Soc., Vol. 101, 1975, pp 13-24 coupled with a clear sky emissivity following Dilley, 1998. `````` Mathias Bavay committed Oct 02, 2012 464 465 466 `````` * If no cloud cover fraction is provided, a parametrization (solar index according to Kasten and Czeplak (1980)) * using iswr_meas and iswr_clear_sky will be used. These parameters can therefore safely * be set to IOUtils::nodata if cloudiness is provided. `````` 467 468 469 470 `````` * @param RH relative humidity (between 0 and 1) * @param TA Air temperature (K) * @param iswr_meas Measured Incoming Short Wave Radiation (W/m^2) * @param iswr_clear_sky Clear Sky Modelled Incoming Short Wave Radiation (W/m^2) `````` Mathias Bavay committed Oct 02, 2012 471 `````` * @param cloudiness Cloud cover fraction (between 0 and 1, optional) `````` 472 473 `````` * @return long wave radiation (W/m^2) or IOUtils::nodata at night time */ `````` Mathias Bavay committed Oct 02, 2012 474 475 476 ``````double Atmosphere::Unsworth_ilwr(const double& RH, const double& TA, const double& iswr_meas, const double& iswr_clear_sky, const double& cloudiness) { double c; `````` 477 `````` if(cloudiness!=IOUtils::nodata) { `````` Mathias Bavay committed Oct 02, 2012 478 479 480 481 482 483 484 485 `````` if(cloudiness<0. || cloudiness>1.) return IOUtils::nodata; c = cloudiness; } else { if(iswr_meas<=0. || iswr_clear_sky<=0.) return IOUtils::nodata; c = Kasten_cloudiness(iswr_meas/iswr_clear_sky); } `````` 486 487 488 `````` const double epsilon_clear = Dilley_emissivity(RH, TA); const double epsilon = (1.-.84*c)*epsilon_clear + .84*c; `````` Mathias Bavay committed Oct 02, 2012 489 `````` `````` 490 491 492 493 494 495 496 `````` return blkBody_Radiation(epsilon, TA); } /** * @brief Evaluate the long wave radiation for clear or cloudy sky. * This uses the formula from Unsworth and Monteith -- "Long-wave radiation at the ground", * Q. J. R. Meteorolo. Soc., Vol. 101, 1975, pp 13-24 coupled with a clear sky emissivity following Dilley, 1998. `````` Mathias Bavay committed Oct 02, 2012 497 498 499 `````` * If no cloud cover fraction is provided, a parametrization (according to Kasten and Czeplak (1980)) * using the current location (lat, lon, altitude) and ISWR will be used. These parameters can therefore safely * be set to IOUtils::nodata if cloudiness is provided. `````` 500 501 502 503 504 505 506 507 `````` * @param lat latitude of the point of observation * @param lon longitude of the point of observation * @param altitude altitude of the point of observation * @param julian julian date at the point of observation * @param TZ time zone at the point of observation * @param RH relative humidity (between 0 and 1) * @param TA Air temperature (K) * @param ISWR Measured Incoming Short Wave Radiation (W/m^2) `````` Mathias Bavay committed Oct 02, 2012 508 `````` * @param cloudiness Cloud cover fraction (between 0 and 1, optional) `````` 509 510 511 512 513 514 515 `````` * @return long wave radiation (W/m^2) or IOUtils::nodata at night time * Please note that this call might NOT be efficient for processing large amounts of points, * since it internally builds complex objects at every call. You might want to copy/paste * its code in order to process data in bulk. */ double Atmosphere::Unsworth_ilwr(const double& lat, const double& lon, const double& altitude, const double& julian, const double& TZ, `````` Mathias Bavay committed Oct 02, 2012 516 `````` const double& RH, const double& TA, const double& ISWR, const double& cloudiness) `````` 517 ``````{ `````` Mathias Bavay committed Oct 02, 2012 518 `````` if(TA==IOUtils::nodata || RH==IOUtils::nodata) { `````` 519 520 521 `````` return IOUtils::nodata; } `````` Mathias Bavay committed Oct 02, 2012 522 523 `````` if(cloudiness==IOUtils::nodata) { if(ISWR==IOUtils::nodata) return IOUtils::nodata; `````` 524 `````` `````` Mathias Bavay committed Oct 02, 2012 525 526 `````` SunObject Sun(lat, lon, altitude, julian, TZ); Sun.calculateRadiation(TA, RH, 0.5); //we force a terrain albedo of 0.5... `````` 527 528 529 `````` double toa, direct, diffuse; Sun.getHorizontalRadiation(toa, direct, diffuse); return Atmosphere::Unsworth_ilwr(RH, TA, ISWR, direct+diffuse); `````` Mathias Bavay committed Oct 02, 2012 530 531 532 `````` } else { return Atmosphere::Unsworth_ilwr(RH, TA, IOUtils::nodata, IOUtils::nodata, cloudiness); } `````` 533 534 ``````} `````` 535 536 537 538 539 540 541 542 543 544 545 546 547 ``````/** * @brief Compute a parametrized Incoming Long Wave Radiation * This uses either Atmosphere::Crawford_ilwr or Atmosphere::Omstedt_ilwr or Atmosphere::Brutsaert_ilwr * depending on which parameters are available. * * @param lat latitude of the point of observation * @param lon longitude of the point of observation * @param altitude altitude of the point of observation * @param julian julian date at the point of observation * @param TZ time zone at the point of observation * @param RH relative humidity (between 0 and 1) * @param TA Air temperature (K) * @param ISWR Measured Incoming Short Wave Radiation (W/m^2) `````` Mathias Bavay committed Oct 02, 2012 548 `````` * @param cloudiness fractional cloud cover (between 0 and 1, optional. If provided, it will have the priority) `````` 549 550 551 552 553 554 555 556 557 `````` * @return long wave radiation (W/m^2) or IOUtils::nodata */ double Atmosphere::ILWR_parametrized(const double& lat, const double& lon, const double& altitude, const double& julian, const double& TZ, const double& RH, const double& TA, const double& ISWR, const double& cloudiness) { const double iswr_thresh = 5.; //any iswr less than this is not considered as valid for Crawford const double ND=IOUtils::nodata; //since we will do lots of comparisons with it... `````` Mathias Bavay committed Oct 02, 2012 558 `````` if(RH!=ND && TA!=ND && cloudiness!=ND) { `````` 559 `````` return Omstedt_ilwr(RH, TA, cloudiness); `````` Mathias Bavay committed Oct 02, 2012 560 561 `````` } if(lat!=ND && lon!=ND && altitude!=ND && julian!=ND && TZ!=ND && RH!=ND && TA!=ND && ISWR!=ND && ISWR>iswr_thresh) { `````` Mathias Bavay committed Jun 14, 2013 562 `````` const double ilwr_p = Unsworth_ilwr(lat, lon, altitude, julian, TZ, RH, TA, ISWR); `````` Fierz committed Oct 26, 2013 563 `````` if(ilwr_p!=ND) return ilwr_p; //it might have been that we could not compute (for low solar angles) `````` Mathias Bavay committed Oct 02, 2012 564 565 `````` } if(RH!=ND && TA!=ND) { `````` 566 567 568 569 570 `````` return Brutsaert_ilwr(RH, TA); } return ND; } `````` 571 `````` `````` Mathias Bavay committed Jan 26, 2011 572 ``````/** `````` Mathias Bavay committed Aug 16, 2011 573 ``````* @brief Convert a relative humidity to a dew point temperature. `````` Mathias Bavay committed Jan 26, 2011 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 ``````* @param RH relative humidity between 0 and 1 * @param TA air temperature (K) * @param force_water if set to true, compute over water. Otherwise, a smooth transition between over ice and over water is computed. * @return dew point temperature (K) */ double Atmosphere::RhtoDewPoint(double RH, double TA, const bool& force_water) { TA = K_TO_C(TA); double Es, E, Tdw, Tdi; //saturation and current water vapor pressure const double Aw = 611.21, Bw = 17.502, Cw = 240.97; //parameters for water const double Ai = 611.15, Bi = 22.452, Ci = 272.55; //parameters for ice const double Tfreeze = 0.; //freezing temperature const double Tnucl = -16.0; //nucleation temperature const double di = 1. / ((TA - Tnucl) * (TA - Tnucl) + 1e-6); //distance to pure ice const double dw = 1. / ((Tfreeze - TA) * (Tfreeze - TA) + 1e-6); //distance to pure water //in order to avoid getting NaN if RH=0 RH += 0.0001; assert(RH>0.); if (TA >= Tfreeze || force_water==true) {//above freezing point, water Es = Aw * exp( (Bw * TA) / (Cw + TA) ); E = RH * Es; Tdw = ( Cw * log(E / Aw) ) / ( Bw - log(E / Aw) ); return C_TO_K(Tdw); } if (TA < Tnucl) { //below nucleation, ice Es = Ai * exp( (Bi * TA) / (Ci + TA) ); E = RH * Es; Tdi = ( Ci * log(E / Ai) ) / ( Bi - log(E / Ai) ); return C_TO_K(Tdi); } //no clear state, we do a smooth interpolation between water and ice Es = Ai * exp( (Bi*TA) / (Ci + TA) ); E = RH * Es; Tdi = ( Ci * log(E / Ai) ) / ( Bi - log(E / Ai) ); Es = Aw * exp( (Bw * TA) / (Cw + TA) ); E = RH * Es; Tdw = ( Cw * log(E / Aw) ) / ( Bw - log(E / Aw) ); return C_TO_K( (di / (di + dw) * Tdi + dw / (di + dw) * Tdw) ); } /** `````` Mathias Bavay committed Aug 16, 2011 619 ``````* @brief Convert a dew point temperature to a relative humidity. `````` Mathias Bavay committed Jan 26, 2011 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 ``````* @param TD dew point temperature (K) * @param TA air temperature (K) * @param force_water if set to true, compute over water. Otherwise, a smooth transition between over ice and over water is computed. * @return relative humidity between 0 and 1 */ double Atmosphere::DewPointtoRh(double TD, double TA, const bool& force_water) { //Convert a dew point temperature into a Relative Humidity //TA, TD are in Kelvins, RH is returned between 0 and 1 TA = K_TO_C(TA); TD = K_TO_C(TD); double Es, E, Rhi, Rhw, Rh; //saturation and current water vapro pressure const double Aw = 611.21, Bw = 17.502, Cw = 240.97; //parameters for water const double Ai = 611.15, Bi = 22.452, Ci = 272.55; //parameters for ice const double Tfreeze = 0.; //freezing temperature const double Tnucl = -16.0; //nucleation temperature const double di = 1. / ((TA - Tnucl) * (TA - Tnucl) + 1e-6); //distance to pure ice const double dw = 1. / ((Tfreeze - TA) * (Tfreeze - TA) + 1e-6); //distance to pure water if (TA >= Tfreeze || force_water==true) { //above freezing point, water Es = Aw * exp( (Bw * TA) / (Cw + TA) ); E = Aw * exp( (Bw * TD) / (Cw + TD) ); Rhw = (E / Es); if (Rhw > 1.) { return 1.; } else { return Rhw; } } if (TA < Tnucl) { //below nucleation, ice Es = Ai * exp( (Bi * TA) / (Ci + TA) ); E = Ai * exp( (Bi * TD) / (Ci + TD) ); Rhi = (E / Es); if (Rhi > 1.) { return 1.; } else { return Rhi; } } //no clear state, we do a smooth interpolation between water and ice Es = Ai * exp( (Bi * TA) / (Ci + TA) ); E = Ai * exp( (Bi * TD) / (Ci + TD) ); Rhi = E / Es; Es = Aw * exp( (Bw * TA) / (Cw + TA) ); E = Aw * exp( (Bw * TD) / (Cw + TD) ); Rhw = E / Es; Rh = (di / (di + dw) * Rhi + dw / (di + dw) * Rhw); if(Rh > 1.) { return 1.; } else { return Rh; } } `````` 679 ``````/** `````` Mathias Bavay committed Aug 16, 2011 680 ``````* @brief Calculate the relative Humidity (RH) from specific humidity. `````` 681 682 ``````* @param altitude altitude over sea level (m) * @param TA air temperature (K) `````` Mathias Bavay committed Aug 16, 2011 683 ``````* @param qi specific humidity `````` 684 685 686 ``````* @return relative humidity between 0 and 1 */ double Atmosphere::specToRelHumidity(const double& altitude, const double& TA, const double& qi) `````` Mathias Bavay committed Aug 16, 2011 687 ``````{ `````` 688 689 690 `````` const double SatVaporDensity = waterVaporDensity(TA, waterSaturationPressure(TA)); const double RH = (qi/(1.-qi))*stdDryAirDensity(altitude, TA)/SatVaporDensity; `````` Mathias Bavay committed Aug 16, 2011 691 692 `````` if(RH>1.) return 1.; else return RH; `````` 693 694 695 ``````} /** `````` Mathias Bavay committed Aug 16, 2011 696 ``````* @brief Calculate the specific Humidity from relative humidity (RH). `````` 697 698 699 700 701 702 703 704 705 706 707 708 709 ``````* @param altitude altitude over sea level (m) * @param TA air temperature (K) * @param RH relative humidity (between 0 and 1) * @return specific humidity */ double Atmosphere::relToSpecHumidity(const double& altitude, const double& TA, const double& RH) { const double dryAir_density = stdDryAirDensity(altitude, TA); const double wetAir_density = RH * waterVaporDensity(TA,waterSaturationPressure(TA)); const double qi = 1./( dryAir_density/wetAir_density+1. ); return qi; } `````` Mathias Bavay committed Jan 26, 2011 710 ``} //namespace``