WSL/SLF GitLab Repository

Atmosphere.cc 37.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/***********************************************************************************/
/*  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 <http://www.gnu.org/licenses/>.
*/

#include <assert.h>
20
#include <sstream>
21
#include <cmath>
22
23
24
#include <meteoio/meteoLaws/Atmosphere.h>
#include <meteoio/meteoLaws/Meteoconst.h>
#include <meteoio/meteoLaws/Sun.h>
25
#include <meteoio/MathOptim.h>
26
27
28
#include <meteoio/IOUtils.h>

namespace mio {
29
30
31
32
/**
 * @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)
33
 * @return black body emissivity (0-1)
34
35
36
 */
double Atmosphere::blkBody_Emissivity(const double& lwr, const double& T) {
	const double T2 = T*T;
37
	const double ea = lwr / (Cst::stefan_boltzmann * (T2*T2));
38
	return std::min(ea, 1.);
39
40
41
42
43
44
}

/**
 * @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)
45
 * @return black body radiation (W/m^2)
46
47
48
49
50
51
 */
double Atmosphere::blkBody_Radiation(const double& ea, const double& T) {
	const double T2 = T*T;
	return ( ea * (Cst::stefan_boltzmann * (T2*T2)) );
}

52
/**
53
54
55
56
57
58
59
* @brief Standard atmospheric pressure as a function of the altitude.
* This is based on the following formula (with h the altitude, P<sub>0</sub> 
* and T<sub>0</sub> the standard sea level pressure and temperature, L the 
* dry adiabatic lapse rate and R<sub>0</sub> the earth's radius):
* \f[
* P = P_0 \cdot {\left( 1 - \frac{L \cdot R_0 \cdot h}{T_0 ( R_0 + h )}  \right)}^{\frac{g}{L R}}
* \f]
60
61
62
63
* @param altitude altitude above sea level (m)
* @return standard pressure (Pa)
*/
double Atmosphere::stdAirPressure(const double& altitude) {
64
65
	const double expo = Cst::gravity / (Cst::mean_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
	const double p = Cst::std_press * pow( 1. - ( (Cst::mean_adiabatique_lapse_rate * Cst::earth_R0 * altitude) / (Cst::std_temp * (Cst::earth_R0 + altitude)) ), expo );
66
67
	return p;
}
68

69
70
71
72
73
74
75
76
77
78
/**
 * @brief Acceleration due to gravity
 * @param altitude altitude above sea level (m)
 * @param latitude latitude in degrees
 * @return acceleration due to gravity (m/s2)
 */
double Atmosphere::gravity(const double& altitude, const double& latitude) {
	const double lat = latitude*Cst::to_rad;
	const double g = 9.780356 * (1. + 0.0052885*Optim::pow2(sin(lat)) - 0.0000059*Optim::pow2(sin(2.*lat))) - 0.003086 * altitude*1e-3;
	return g;
79
80
}

81
82
83
84
85
86
87
88
89
90
91
92
/**
* @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
93
* @param Temperature air temperature (K)
94
95
96
97
98
99
100
* @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);
}

101
/**
102
* @brief Standard atmosphere wet bulb temperature.
103
104
105
106
107
108
109
110
111
112
113
114
* 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);
115
	const double Vp = vaporSaturationPressure(T);
116
117
118
119

	return ( T - (RH*Vp - Vp) * mixing_ratio * L / p / Cst::specific_heat_air );
}

120
121
122
/**
* @brief Black Globe Temperature.
* This is an estimation of the black globe temperature based on physical modeling as in
123
124
* V. E. Dimiceli, S. F. Piltz and S. A. Amburn, <i>"Estimation of Black Globe Temperature for Calculation of the 
* Wet Bulb Globe Temperature Index"</i> in World Congress on Engineering and Computer Science, <b>2</b>, 2011.
125
126
127
128
129
130
131
132
133
134
135
* @param TA air temperature (K)
* @param RH relative humidity (between 0 and 1)
* @param VW wind velocity (m/s)
* @param iswr_dir direct solar SW radiation (W/m^2)
* @param iswr_diff diffuse solar SW radiation (W/m^2)
* @param cos_Z cosinus of the solar zenith angle
* @return black globe temperature (K)
*/
double Atmosphere::blackGlobeTemperature(const double& TA, const double& RH, const double& VW, const double& iswr_dir, const double& iswr_diff, const double& cos_Z)
{
	const double S = iswr_dir + iswr_diff;
136
	//const double a=1, b=1, c=1; // get real values!
137
138
	//const double h = a * pow(S, b) * pow(cos_Z, c);
	const double h = 0.315; //personnal communication from S. Amburn
139
	const double emissivity = 0.575 * pow(RH*vaporSaturationPressure(TA), 1./7.);
140
141
	const double B = S * (iswr_dir/(4.*Cst::stefan_boltzmann*cos_Z) + 1.2/Cst::stefan_boltzmann*iswr_diff) + emissivity*Optim::pow4(TA);
	const double C = h * pow(VW, 0.58) / 5.3865e-8;
142

143
	const double Tg = (B + C*TA + 7680000.) / (C + 256000.);
144
145
146
	return Tg;
}

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
/**
* @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:
* <center><table border="0">
* <tr><th>land class</th><th>roughness length</th></tr>
* <tr><td>sea</td><td>0.0002</td></tr>
* <tr><td>smooth</td><td>0.005</td></tr>
* <tr><td>open</td><td>0.03</td></tr>
* <tr><td>roughly open</td><td>0.10</td></tr>
* <tr><td>rough</td><td>0.25</td></tr>
* <tr><td>very rough</td><td>0.5</td></tr>
* <tr><td>closed</td><td>1.0</td></tr>
* <tr><td>chaotic</td><td>over 2.0</td></tr>
* </table></center>
*
* @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) );
}

175
176
177
178
179
180
181
182
183
184
185
/**
* @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)
{
186
	if (TA>(273.15+10.) || VW<5.) return TA; //not applicable in this case
187

188
	const double t = IOUtils::K_TO_C(TA); //in Celsius
189
190
191
192
	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;
193
	return IOUtils::C_TO_K(WCT);
194
195
196
197
198
199
200
201
202
203
204
205
206
}

/**
* @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)
{
207
	if (TA<(273.15+26.7) || RH<0.4) return TA; //not applicable in this case
208

209
	const double t = IOUtils::K_TO_C(TA); //in Celsius
210
211
212
213
214
215
216
217
218
	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;

219
	return IOUtils::C_TO_K(HI);
220
221
}

222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
/**
* @brief Wet Bulb Globe Temperature index.
* This is an index aiming at expressing the human-perceived air temperature due to humidity, wind and radiation.
* This is the foundation of ISO 7243 and is widely used for physical activity safety evaluation (for example for physical training).
* @param TA air temperature (K)
* @param RH relative humidity (between 0 and 1)
* @param VW wind velocity (m/s)
* @param iswr_dir direct solar SW radiation (W/m^2)
* @param iswr_diff diffuse solar SW radiation (W/m^2)
* @param cos_Z cosinus of the solar zenith angle
* @param altitude altitude of the point where to perform the calculation
* @return Heat index (K)
*/
double Atmosphere::WBGT_index(const double& TA, const double& RH, const double& VW, const double& iswr_dir, const double& iswr_diff, const double& cos_Z, const double& altitude)
{
	const double NWB = wetBulbTemperature(TA, RH, altitude);
	const double GT = blackGlobeTemperature(TA, RH, VW, iswr_dir, iswr_diff, cos_Z);
	const double DB = TA;
240

241
242
243
	return 0.7*NWB + 0.2*GT + 0.1*DB;
}

244
/**
245
246
247
* @brief Standard water vapor saturation pressure.
* See Murray, F. W., <i>"On the computation of saturation vapor pressure"</i>, 1966, J. Appl. Meteor., <b>6</b>, 203–204,
* doi: 10.1175/1520-0450(1967)006<0203:OTCOSV>2.0.CO;2.
248
249
250
* @param T air temperature (K)
* @return standard water vapor saturation pressure (Pa)
*/
251
double Atmosphere::vaporSaturationPressure(const double& T) {
252
253
	double c2, c3; // varying constants

254
	if ( T < Cst::t_water_triple_pt ) { // for a flat ice surface
255
256
257
258
259
260
		c2 = 21.88;
		c3 = 7.66;
	} else { // for a flat water surface
		c2 = 17.27;
		c3 = 35.86;
	}
261
	const double exp_p_sat = c2 *  (T - Cst::t_water_triple_pt) / (T - c3); //exponent
262
263
264
265

	return( Cst::p_water_triple_pt * exp( exp_p_sat ) );
}

266
267
268
269
270
271
272
273
274
275
276
277
278
/**
* @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);
}

279
/**
280
 * @brief Evaluate the atmosphere emissivity for clear sky.
281
282
 * This uses the formula from Brutsaert -- "On a Derivable
 * Formula for Long-Wave Radiation From Clear Skies", Journal of Water Resources
283
 * Research, <b>11</b>, No. 5, October 1975, pp 742-744.
284
 * Alternative: Satterlund (1979): Water Resources Research, 15, 1649-1650.
285
 * @param RH relative humidity (between 0 and 1)
286
287
288
 * @param TA Air temperature (K)
 * @return clear sky emissivity
 */
289
double Atmosphere::Brutsaert_emissivity(const double& RH, const double& TA) {
290
	const double e0 = RH * vaporSaturationPressure(TA); //water vapor pressure
291
292
	const double e0_mBar = 0.01 * e0;
	const double exponent = 1./7.;
293
	const double ea = 1.24 * pow( (e0_mBar / TA), exponent);
294
	return std::min(ea, 1.);
295
296
297
}

/**
298
 * @brief Evaluate the long wave radiation for clear sky.
299
300
 * This uses the formula from Brutsaert -- "On a Derivable
 * Formula for Long-Wave Radiation From Clear Skies", Journal of Water Resources
301
 * Research, <b>11</b>, No. 5, October 1975, pp 742-744.
302
303
304
305
306
307
 * 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) {
308
	const double ea = Brutsaert_emissivity(RH, TA);
309
	return blkBody_Radiation(ea, TA);
310
311
}

312
313
314
315
/**
 * @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",
316
 * Q. J. R. Meteorolo. Soc., <b>124</b>, 1998, pp 1391-1401. The long wave is computed
317
318
319
 * 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)
320
 * @return clear sky emissivity
321
322
323
324
*/
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);
325
	const double ea = ilwr_dilley/ilwr_blkbody;
326
	return std::min(ea, 1.);
327
328
329
330
331
332
}

/**
 * @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",
333
 * Q. J. R. Meteorolo. Soc., <b>124</b>, 1998, pp 1391-1401.
334
335
336
337
338
 * @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) {
339
	const double e0 = RH * vaporSaturationPressure(TA) * 0.001; //water vapor pressure, kPa
340
341
342
343
344
345
346
347
348
349
	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
350
 * downward clear-sky radiation at the surface", Q. J. R. Meteorolo. Soc., <b>122</b>, 1996, pp 1127-1151.
351
352
 * @param RH relative humidity (between 0 and 1)
 * @param TA near surface air temperature (K)
353
 * @return clear sky emissivity
354
355
*/
double Atmosphere::Prata_emissivity(const double& RH, const double& TA) {
356
	const double e0 = RH * vaporSaturationPressure(TA) * 0.001; //water vapor pressure, kPa
357
	const double w = 4650.*e0/TA; //precipitable water, Prata 1996
358
	const double ea = 1. - (1.+w)*exp( -sqrt(1.2+3.*w) );
359
	return std::min(ea, 1.);
360
361
362
363
364
}

/**
 * @brief Evaluate the long wave radiation for clear sky.
 * This uses the formula from Prata -- "A new long-wave formula for estimating
365
 * downward clear-sky radiation at the surface", Q. J. R. Meteorolo. Soc., <b>122</b>, 1996, pp 1127-1151.
366
367
368
369
370
371
372
 * @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);
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
}

/**
 * @brief Evaluate the atmosphere emissivity for clear sky.
 * This uses the formula from Clark & Allen -- "The estimation of atmospheric radiation for clear and
 * cloudy skies", Proceedings of the second national passive solar conference, <b>2</b>, 1978, p 676.
 * @param RH relative humidity (between 0 and 1)
 * @param TA near surface air temperature (K)
 * @return clear sky emissivity
*/
double Atmosphere::Clark_emissivity(const double& RH, const double& TA) {
	const double Tdp = RhtoDewPoint(RH, TA, false);
	const double ea = 0.787 + 0.0028 * (Tdp-Cst::t_water_triple_pt);
	return ea;
}

/**
 * @brief Evaluate the long wave radiation for clear sky.
 * This uses the formula from Clark & Allen -- "The estimation of atmospheric radiation for clear and
 * cloudy skies", Proceedings of the second national passive solar conference, <b>2</b>, 1978, p 676.
 * @param RH relative humidity (between 0 and 1)
 * @param TA near surface air temperature (K)
 * @return long wave radiation (W/m^2)
*/
double Atmosphere::Clark_ilwr(const double& RH, const double& TA) {
	const double epsilon = Clark_emissivity(RH, TA);
	return blkBody_Radiation(epsilon, TA);
}

/**
 * @brief Evaluate the atmosphere emissivity for clear sky.
 * This uses the formula from Tang, Etzion and Meir -- "Estimates of clear night sky emissivity in the
 * Negev Highlands, Israel", Energy Conversion and Management, <b>45.11</b>, 2004, pp 1831-1843.
 * @param RH relative humidity (between 0 and 1)
 * @param TA near surface air temperature (K)
 * @return clear sky emissivity
*/
double Atmosphere::Tang_emissivity(const double& RH, const double& TA) {
	const double Tdp = RhtoDewPoint(RH, TA, false);
	const double ea = 0.754 + 0.0044 * (Tdp-Cst::t_water_triple_pt);
413
	return std::min(ea, 1.);
414
415
416
417
418
419
420
421
422
423
424
425
426
}

/**
 * @brief Evaluate the long wave radiation for clear sky.
 * This uses the formula from Tang, Etzion and Meir -- "Estimates of clear night sky emissivity in the
 * Negev Highlands, Israel", Energy Conversion and Management, <b>45.11</b>, 2004, pp 1831-1843.
 * @param RH relative humidity (between 0 and 1)
 * @param TA near surface air temperature (K)
 * @return long wave radiation (W/m^2)
*/
double Atmosphere::Tang_ilwr(const double& RH, const double& TA) {
	const double epsilon = Tang_emissivity(RH, TA);
	return blkBody_Radiation(epsilon, TA);
427
428
429
430
431
432
433
434
435
436
437
}

/**
 * @brief Evaluate the atmosphere emissivity for clear sky.
 * This uses the formula from Idso -- "A set of equations for full spectrum and 8 to 14 um and
 * 10.5 to 12.5 um thermal radiation from cloudless skies", Water Resources Research, <b>17</b>, 1981, pp 295-304.
 * @param RH relative humidity (between 0 and 1)
 * @param TA near surface air temperature (K)
 * @return clear sky emissivity
*/
double Atmosphere::Idso_emissivity(const double& RH, const double& TA) {
438
	const double e0 = RH * vaporSaturationPressure(TA) * 0.0001; //water vapor pressure, mbar
439
	const double ea = 0.70 + 5.95e-5 * e0 * exp(1500./TA);
440
	return std::min(ea, 1.);
441
442
443
444
445
446
447
448
449
450
451
452
453
}

/**
 * @brief Evaluate the long wave radiation for clear sky.
 * This uses the formula from Idso -- "A set of equations for full spectrum and 8 to 14 um and
 * 10.5 to 12.5 um thermal radiation from cloudless skies", Water Resources Research, <b>17</b>, 1981, pp 295-304.
 * @param RH relative humidity (between 0 and 1)
 * @param TA near surface air temperature (K)
 * @return long wave radiation (W/m^2)
*/
double Atmosphere::Idso_ilwr(const double& RH, const double& TA) {
	const double epsilon = Idso_emissivity(RH, TA);
	return blkBody_Radiation(epsilon, TA);
454
455
}

456
457
458
/**
* @brief Evaluate the atmosphere emissivity from the water vapor pressure and cloudiness.
* This is according to A. Omstedt, <i>"A coupled one-dimensional sea ice-ocean model applied to a semi-enclosed basin"</i>,
459
* Tellus, <b>42 A</b>, 568-582, 1990, DOI:10.1034/j.1600-0870.1990.t01-3-00007.x
460
461
462
463
464
465
* @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 emissivity (between 0 and 1)
*/
double Atmosphere::Omstedt_emissivity(const double& RH, const double& TA, const double& cloudiness) {
466
	const double e0 = RH * vaporSaturationPressure(TA); //water vapor pressure
467
468
469
470
471
472
	const double eps_w = 0.97;
	const double a1 = 0.68;
	const double a2 = 0.0036;
	const double a3 = 0.18;

	const double ea = (eps_w * (a1 + a2 * sqrt(e0)) * (1. + a3 * cloudiness * cloudiness)); //emissivity
473
	return std::min(ea, 1.);
474
475
476
477
478
}

/**
* @brief Evaluate the long wave radiation from RH, TA and cloudiness.
* This is according to A. Omstedt, <i>"A coupled one-dimensional sea ice-ocean model applied to a semi-enclosed basin"</i>,
479
* Tellus, <b>42 A</b>, 568-582, 1990, DOI:10.1034/j.1600-0870.1990.t01-3-00007.x
480
481
482
483
484
485
486
487
488
489
* @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) {
	const double ea = Omstedt_emissivity(RH, TA, cloudiness);
	return blkBody_Radiation(ea, TA);
}

490
491
492
493
494
495
496
497
498
499
/**
* @brief Evaluate the atmosphere emissivity from the water vapor pressure and cloudiness.
* This is according to Konzelmann, Thomas, et al. <i>"Parameterization of global and longwave incoming radiation
* for the Greenland Ice Sheet."</i> Global and Planetary change <b>9.1</b> (1994): 143-164.
* @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 emissivity (between 0 and 1)
*/
double Atmosphere::Konzelmann_emissivity(const double& RH, const double& TA, const double& cloudiness) {
500
	const double ea = RH * vaporSaturationPressure(TA); //screen-level water vapor pressure
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
	const double exponent = 1./8.;

	const double epsilon_cs = 0.23 + 0.484*pow( ea/TA, exponent ); //clear sky emissivity
	const double epsilon_oc = 0.952; //fully overcast sky emissivity

	const double weight = Optim::pow4(cloudiness); //weight for the weighted average between clear sky and overcast

	const double epsilon_cloudy = (1. - weight)*epsilon_cs + weight*epsilon_oc;
	return epsilon_cloudy;
}

/**
* @brief Evaluate the long wave radiation from RH, TA and cloudiness.
* This is according to Konzelmann, Thomas, et al. <i>"Parameterization of global and longwave incoming radiation
* for the Greenland Ice Sheet."</i> Global and Planetary change <b>9.1</b> (1994): 143-164.
* @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::Konzelmann_ilwr(const double& RH, const double& TA, const double& cloudiness) {
	const double ea = Konzelmann_emissivity(RH, TA, cloudiness);
	return blkBody_Radiation(ea, TA);
}

526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
/**
 * @brief Evaluate the solar clearness index for a given cloudiness. 
 * This uses the formula from Kasten and Czeplak -- <i>"Solar and terrestrial radiation
 * dependent on the amount and type of cloud"</i>, Sol. Energy, <b>24</b>, 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 cloudiness in okta, between 0 and 1
 * @return solar clearness index
*/
double Atmosphere::Kasten_clearness(const double& cloudiness) {
	const double b1 = 0.75, b2 = 3.4;
	if (cloudiness<0. || cloudiness>1.) {
		std::ostringstream ss;
		ss << "Invalid cloudiness value: " << cloudiness << " (it should be between 0 and 1)";
		throw InvalidArgumentException(ss.str(), AT);
	}
	const double clearness = 1. - b1*pow(cloudiness, b2);
	return clearness;
}

546
547
548
549
550
551
552
/**
 * @brief Evaluate the cloudiness from a given solar index.
 * This uses the formula from Kasten and Czeplak -- <i>"Solar and terrestrial radiation
 * dependent on the amount and type of cloud"</i>, Sol. Energy, <b>24</b>, 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
553
 * @return cloudiness (in okta, between 0 and 1)
554
555
556
557
*/
double Atmosphere::Kasten_cloudiness(const double& solarIndex) {
	const double b1 = 0.75, b2 = 3.4;

558
	if (solarIndex>1.) return 0.;
559
	const double cloudiness = pow((1.-solarIndex)/b1, 1./b2);
560
	return std::min(cloudiness, 1.);
561
562
}

563
564
565
566
567
/**
 * @brief Evaluate the long wave radiation for clear or cloudy sky.
 * This uses the formula from Crawford and Duchon -- <i>"An Improved Parametrization
 * for Estimating Effective Atmospheric Emissivity for Use in Calculating Daytime
 * Downwelling Longwave Radiation"</i>, Journal of Applied Meteorology,
568
569
570
 * <b>38</b>, 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.
571
572
573
574
575
 * @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)
576
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
577
578
 * @return long wave radiation (W/m^2) or IOUtils::nodata at night time
*/
579
580
581
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;
582
583
	if (cloudiness==IOUtils::nodata) {
		if (iswr_meas<=0. || iswr_clear_sky<=0.)
584
585
			return IOUtils::nodata;
		clf = 1. - iswr_meas/iswr_clear_sky;  //cloud fraction estimate
586
		if (clf<0.) clf=0.;
587
	} else {
588
		if (cloudiness<0. || cloudiness>1.)
589
590
591
			return IOUtils::nodata;
		clf = cloudiness;
	}
592

593
	const double e = RH * vaporSaturationPressure(TA); //near surface water vapor pressure
594
595
596
	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.);
597
	return blkBody_Radiation(epsilon, TA);
598
599
600
601
602
603
604
}

/**
 * @brief Evaluate the long wave radiation for clear or cloudy sky.
 * This uses the formula from Crawford and Duchon -- <i>"An Improved Parametrization
 * for Estimating Effective Atmospheric Emissivity for Use in Calculating Daytime
 * Downwelling Longwave Radiation"</i>, Journal of Applied Meteorology,
605
606
607
 * <b>38</b>, 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.
608
609
610
611
612
613
614
615
 * @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)
616
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
617
618
619
620
621
622
623
 * @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,
624
                                 const double& RH, const double& TA, const double& ISWR, const double& cloudiness)
625
{
626
	if (TA==IOUtils::nodata || RH==IOUtils::nodata) {
627
628
629
630
631
632
633
		return IOUtils::nodata;
	}

	Date date(julian, TZ, 0.);
	int year, month, day;
	date.getDate(year, month, day);

634
635
	if (cloudiness==IOUtils::nodata) {
		if (ISWR==IOUtils::nodata) return IOUtils::nodata;
636
637
638

		SunObject Sun(lat, lon, altitude, julian, TZ);
		Sun.calculateRadiation(TA, RH, 0.5); //we force a terrain albedo of 0.5...
639
640
641
		double toa, direct, diffuse;
		Sun.getHorizontalRadiation(toa, direct, diffuse);
		return Atmosphere::Crawford_ilwr(RH, TA, ISWR, direct+diffuse, static_cast<unsigned char>(month));
642
643
644
	} else {
		return Atmosphere::Crawford_ilwr(RH, TA, IOUtils::nodata, IOUtils::nodata, static_cast<unsigned char>(month), cloudiness);
	}
645
646
}

647
648
649
650
/**
 * @brief Evaluate the long wave radiation for clear or cloudy sky.
 * This uses the formula from Unsworth and Monteith -- <i>"Long-wave radiation at the ground"</i>,
 * Q. J. R. Meteorolo. Soc., Vol. 101, 1975, pp 13-24 coupled with a clear sky emissivity following Dilley, 1998.
651
652
653
 * 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.
654
655
656
657
 * @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)
658
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
659
660
 * @return long wave radiation (W/m^2) or IOUtils::nodata at night time
*/
661
662
663
double Atmosphere::Unsworth_ilwr(const double& RH, const double& TA, const double& iswr_meas, const double& iswr_clear_sky, const double& cloudiness)
{
	double c;
664
665
	if (cloudiness!=IOUtils::nodata) {
		if (cloudiness<0. || cloudiness>1.)
666
667
668
			return IOUtils::nodata;
		c = cloudiness;
	} else {
669
		if (iswr_meas<=0. || iswr_clear_sky<=0.)
670
671
672
			return IOUtils::nodata;
		c = Kasten_cloudiness(iswr_meas/iswr_clear_sky);
	}
673
674
675

	const double epsilon_clear = Dilley_emissivity(RH, TA);
	const double epsilon = (1.-.84*c)*epsilon_clear + .84*c;
676

677
678
679
680
681
682
683
	return blkBody_Radiation(epsilon, TA);
}

/**
 * @brief Evaluate the long wave radiation for clear or cloudy sky.
 * This uses the formula from Unsworth and Monteith -- <i>"Long-wave radiation at the ground"</i>,
 * Q. J. R. Meteorolo. Soc., Vol. 101, 1975, pp 13-24 coupled with a clear sky emissivity following Dilley, 1998.
684
685
686
 * 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.
687
688
689
690
691
692
693
694
 * @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)
695
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
696
697
698
699
700
701
702
 * @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,
703
                                 const double& RH, const double& TA, const double& ISWR, const double& cloudiness)
704
{
705
	if (TA==IOUtils::nodata || RH==IOUtils::nodata) {
706
707
708
		return IOUtils::nodata;
	}

709
710
	if (cloudiness==IOUtils::nodata) {
		if (ISWR==IOUtils::nodata) return IOUtils::nodata;
711

712
713
		SunObject Sun(lat, lon, altitude, julian, TZ);
		Sun.calculateRadiation(TA, RH, 0.5); //we force a terrain albedo of 0.5...
714
715
716
		double toa, direct, diffuse;
		Sun.getHorizontalRadiation(toa, direct, diffuse);
		return Atmosphere::Unsworth_ilwr(RH, TA, ISWR, direct+diffuse);
717
718
719
	} else {
		return Atmosphere::Unsworth_ilwr(RH, TA, IOUtils::nodata, IOUtils::nodata, cloudiness);
	}
720
721
}

722
723
724
725
726
727
728
729
730
731
732
733
734
/**
 * @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)
735
 * @param cloudiness fractional cloud cover (between 0 and 1, optional. If provided, it will have the priority)
736
737
738
739
740
741
742
743
744
 * @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...

745
	if (RH!=ND && TA!=ND && cloudiness!=ND) {
746
		return Omstedt_ilwr(RH, TA, cloudiness);
747
	}
748
	if (lat!=ND && lon!=ND && altitude!=ND && julian!=ND && TZ!=ND && RH!=ND && TA!=ND && ISWR!=ND && ISWR>iswr_thresh) {
749
		const double ilwr_p = Unsworth_ilwr(lat, lon, altitude, julian, TZ, RH, TA, ISWR);
750
		if (ilwr_p!=ND) return ilwr_p; //it might have been that we could not compute (for low solar angles)
751
	}
752
	if (RH!=ND && TA!=ND) {
753
754
755
756
757
		return Brutsaert_ilwr(RH, TA);
	}

	return ND;
}
758

759
/**
760
* @brief Convert a relative humidity to a dew point temperature.
761
762
763
764
765
766
767
* @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)
{
768
	TA = IOUtils::K_TO_C(TA);
769
770
771
772
773
774
775
776
777
778
779
780
781
	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

	//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) );
782
		return IOUtils::C_TO_K(Tdw);
783
784
785
786
787
	}
	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) );
788
		return IOUtils::C_TO_K(Tdi);
789
790
791
792
793
794
795
796
797
798
799
	}

	//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) );

800
801
	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
802
	return IOUtils::C_TO_K( (di / (di + dw) * Tdi + dw / (di + dw) * Tdw) );
803
804
805
}

/**
806
* @brief Convert a dew point temperature to a relative humidity.
807
808
809
810
811
812
813
814
815
* @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
816
817
	TA = IOUtils::K_TO_C(TA);
	TD = IOUtils::K_TO_C(TD);
818
819
820
821
	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
822
	double Es, E, Rhi, Rhw;                         //saturation and current water vapro pressure
823
824
825
826
827
828

	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);
829
		return std::min(Rhw, 1.);
830
831
832
833
834
835
	}
	if (TA < Tnucl) {
		//below nucleation, ice
		Es = Ai * exp( (Bi * TA) / (Ci + TA) );
		E  = Ai * exp( (Bi * TD) / (Ci + TD) );
		Rhi = (E / Es);
836
		return std::min(Rhi, 1.);
837
838
839
840
841
842
843
844
845
846
847
	}

	//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;

848
849
850
851
	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
	const double Rh = (di / (di + dw) * Rhi + dw / (di + dw) * Rhw);
	return std::min(Rh, 1.);
852
853
}

854
/**
855
* @brief Calculate the relative Humidity (RH) from specific humidity.
856
857
* @param altitude altitude over sea level (m)
* @param TA air temperature (K)
858
* @param qi specific humidity
859
860
861
* @return relative humidity between 0 and 1
*/
double Atmosphere::specToRelHumidity(const double& altitude, const double& TA, const double& qi)
862
{
863
	const double SatVaporDensity = waterVaporDensity(TA, vaporSaturationPressure(TA));
864
865
	const double dryAir_density = stdDryAirDensity(altitude, TA);
	const double RH = qi/(1.-qi) * dryAir_density/SatVaporDensity;
866

867
	return std::min(RH, 1.);
868
869
870
}

/**
871
* @brief Calculate the specific Humidity from relative humidity (RH).
872
873
874
875
876
877
878
879
* @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);
880
	const double SatVaporDensity = waterVaporDensity(TA, vaporSaturationPressure(TA));
881
882
883
	const double qi_inv = dryAir_density/(RH*SatVaporDensity) + 1.;

	return 1./qi_inv;
884
885
}

886
} //namespace