WSL/SLF GitLab Repository

Atmosphere.cc 37.3 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
115
116
117
118
119
* 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 );
}

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
140
141
	const double emissivity = 0.575 * pow(RH*waterSaturationPressure(TA), 1./7.);
	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
248
249
250
251
/**
* @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

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

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

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

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

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

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

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

/**
 * @brief Evaluate the long wave radiation for clear sky.
 * This uses the formula from Prata -- "A new long-wave formula for estimating
363
 * downward clear-sky radiation at the surface", Q. J. R. Meteorolo. Soc., <b>122</b>, 1996, pp 1127-1151.
364
365
366
367
368
369
370
 * @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);
371
372
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
}

/**
 * @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);
411
	return std::min(ea, 1.);
412
413
414
415
416
417
418
419
420
421
422
423
424
}

/**
 * @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);
425
426
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) {
	const double e0 = RH * waterSaturationPressure(TA) * 0.0001; //water vapor pressure, mbar
	const double ea = 0.70 + 5.95e-5 * e0 * exp(1500./TA);
438
	return std::min(ea, 1.);
439
440
441
442
443
444
445
446
447
448
449
450
451
}

/**
 * @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);
452
453
}

454
455
456
/**
* @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>,
457
* Tellus, <b>42 A</b>, 568-582, 1990, DOI:10.1034/j.1600-0870.1990.t01-3-00007.x
458
459
460
461
462
463
464
465
466
467
468
469
470
* @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) {
	const double e0 = RH * waterSaturationPressure(TA); //water vapor pressure
	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
471
	return std::min(ea, 1.);
472
473
474
475
476
}

/**
* @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>,
477
* Tellus, <b>42 A</b>, 568-582, 1990, DOI:10.1034/j.1600-0870.1990.t01-3-00007.x
478
479
480
481
482
483
484
485
486
487
* @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);
}

488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
/**
* @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) {
	const double ea = RH * waterSaturationPressure(TA); //screen-level water vapor pressure
	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);
}

524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
/**
 * @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;
}

544
545
546
547
548
549
550
/**
 * @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
551
 * @return cloudiness (in okta, between 0 and 1)
552
553
554
555
*/
double Atmosphere::Kasten_cloudiness(const double& solarIndex) {
	const double b1 = 0.75, b2 = 3.4;

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

561
562
563
564
565
/**
 * @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,
566
567
568
 * <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.
569
570
571
572
573
 * @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)
574
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
575
576
 * @return long wave radiation (W/m^2) or IOUtils::nodata at night time
*/
577
578
579
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;
580
581
	if (cloudiness==IOUtils::nodata) {
		if (iswr_meas<=0. || iswr_clear_sky<=0.)
582
583
			return IOUtils::nodata;
		clf = 1. - iswr_meas/iswr_clear_sky;  //cloud fraction estimate
584
		if (clf<0.) clf=0.;
585
	} else {
586
		if (cloudiness<0. || cloudiness>1.)
587
588
589
			return IOUtils::nodata;
		clf = cloudiness;
	}
590
591
592
593
594

	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.);
595
	return blkBody_Radiation(epsilon, TA);
596
597
598
599
600
601
602
}

/**
 * @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,
603
604
605
 * <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.
606
607
608
609
610
611
612
613
 * @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)
614
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
615
616
617
618
619
620
621
 * @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,
622
                                 const double& RH, const double& TA, const double& ISWR, const double& cloudiness)
623
{
624
	if (TA==IOUtils::nodata || RH==IOUtils::nodata) {
625
626
627
628
629
630
631
		return IOUtils::nodata;
	}

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

632
633
	if (cloudiness==IOUtils::nodata) {
		if (ISWR==IOUtils::nodata) return IOUtils::nodata;
634
635
636

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

645
646
647
648
/**
 * @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.
649
650
651
 * 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.
652
653
654
655
 * @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)
656
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
657
658
 * @return long wave radiation (W/m^2) or IOUtils::nodata at night time
*/
659
660
661
double Atmosphere::Unsworth_ilwr(const double& RH, const double& TA, const double& iswr_meas, const double& iswr_clear_sky, const double& cloudiness)
{
	double c;
662
663
	if (cloudiness!=IOUtils::nodata) {
		if (cloudiness<0. || cloudiness>1.)
664
665
666
			return IOUtils::nodata;
		c = cloudiness;
	} else {
667
		if (iswr_meas<=0. || iswr_clear_sky<=0.)
668
669
670
			return IOUtils::nodata;
		c = Kasten_cloudiness(iswr_meas/iswr_clear_sky);
	}
671
672
673

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

675
676
677
678
679
680
681
	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.
682
683
684
 * 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.
685
686
687
688
689
690
691
692
 * @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)
693
 * @param cloudiness Cloud cover fraction (between 0 and 1, optional)
694
695
696
697
698
699
700
 * @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,
701
                                 const double& RH, const double& TA, const double& ISWR, const double& cloudiness)
702
{
703
	if (TA==IOUtils::nodata || RH==IOUtils::nodata) {
704
705
706
		return IOUtils::nodata;
	}

707
708
	if (cloudiness==IOUtils::nodata) {
		if (ISWR==IOUtils::nodata) return IOUtils::nodata;
709

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

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

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

	return ND;
}
756

757
/**
758
* @brief Convert a relative humidity to a dew point temperature.
759
760
761
762
763
764
765
* @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)
{
766
	TA = IOUtils::K_TO_C(TA);
767
768
769
770
771
772
773
774
775
776
777
778
779
	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) );
780
		return IOUtils::C_TO_K(Tdw);
781
782
783
784
785
	}
	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) );
786
		return IOUtils::C_TO_K(Tdi);
787
788
789
790
791
792
793
794
795
796
797
	}

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

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

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

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

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

846
847
848
849
	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.);
850
851
}

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

865
	return std::min(RH, 1.);
866
867
868
}

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

	return 1./qi_inv;
882
883
}

884
} //namespace