|
1
|
|
|
# Copyright (c) 2008-2015 MetPy Developers. |
|
2
|
|
|
# Distributed under the terms of the BSD 3-Clause License. |
|
3
|
|
|
# SPDX-License-Identifier: BSD-3-Clause |
|
4
|
|
|
"""Contains a collection of thermodynamic calculations.""" |
|
5
|
|
|
|
|
6
|
|
|
from __future__ import division |
|
7
|
|
|
|
|
8
|
|
|
import numpy as np |
|
9
|
|
|
import scipy.integrate as si |
|
10
|
|
|
import scipy.optimize as so |
|
11
|
|
|
|
|
12
|
|
|
from .tools import find_intersections |
|
13
|
|
|
from ..constants import Cp_d, epsilon, kappa, Lv, P0, Rd |
|
14
|
|
|
from ..package_tools import Exporter |
|
15
|
|
|
from ..units import atleast_1d, check_units, concatenate, units |
|
16
|
|
|
|
|
17
|
|
|
exporter = Exporter(globals()) |
|
18
|
|
|
|
|
19
|
|
|
sat_pressure_0c = 6.112 * units.millibar |
|
20
|
|
|
|
|
21
|
|
|
|
|
22
|
|
|
@exporter.export |
|
23
|
|
|
@check_units('[pressure]', '[temperature]') |
|
24
|
|
|
def potential_temperature(pressure, temperature): |
|
25
|
|
|
r"""Calculate the potential temperature. |
|
26
|
|
|
|
|
27
|
|
|
Uses the Poisson equation to calculation the potential temperature |
|
28
|
|
|
given `pressure` and `temperature`. |
|
29
|
|
|
|
|
30
|
|
|
Parameters |
|
31
|
|
|
---------- |
|
32
|
|
|
pressure : `pint.Quantity` |
|
33
|
|
|
The total atmospheric pressure |
|
34
|
|
|
temperature : `pint.Quantity` |
|
35
|
|
|
The temperature |
|
36
|
|
|
|
|
37
|
|
|
Returns |
|
38
|
|
|
------- |
|
39
|
|
|
`pint.Quantity` |
|
40
|
|
|
The potential temperature corresponding to the the temperature and |
|
41
|
|
|
pressure. |
|
42
|
|
|
|
|
43
|
|
|
See Also |
|
44
|
|
|
-------- |
|
45
|
|
|
dry_lapse |
|
46
|
|
|
|
|
47
|
|
|
Notes |
|
48
|
|
|
----- |
|
49
|
|
|
Formula: |
|
50
|
|
|
|
|
51
|
|
|
.. math:: \Theta = T (P_0 / P)^\kappa |
|
52
|
|
|
|
|
53
|
|
|
Examples |
|
54
|
|
|
-------- |
|
55
|
|
|
>>> from metpy.units import units |
|
56
|
|
|
>>> metpy.calc.potential_temperature(800. * units.mbar, 273. * units.kelvin) |
|
57
|
|
|
<Quantity(290.96653180346203, 'kelvin')> |
|
58
|
|
|
|
|
59
|
|
|
""" |
|
60
|
|
|
return temperature * (P0 / pressure).to('dimensionless')**kappa |
|
61
|
|
|
|
|
62
|
|
|
|
|
63
|
|
|
@exporter.export |
|
64
|
|
|
@check_units('[pressure]', '[temperature]') |
|
65
|
|
|
def dry_lapse(pressure, temperature): |
|
66
|
|
|
r"""Calculate the temperature at a level assuming only dry processes. |
|
67
|
|
|
|
|
68
|
|
|
This function lifts a parcel starting at `temperature`, conserving |
|
69
|
|
|
potential temperature. The starting pressure should be the first item in |
|
70
|
|
|
the `pressure` array. |
|
71
|
|
|
|
|
72
|
|
|
Parameters |
|
73
|
|
|
---------- |
|
74
|
|
|
pressure : `pint.Quantity` |
|
75
|
|
|
The atmospheric pressure level(s) of interest |
|
76
|
|
|
temperature : `pint.Quantity` |
|
77
|
|
|
The starting temperature |
|
78
|
|
|
|
|
79
|
|
|
Returns |
|
80
|
|
|
------- |
|
81
|
|
|
`pint.Quantity` |
|
82
|
|
|
The resulting parcel temperature at levels given by `pressure` |
|
83
|
|
|
|
|
84
|
|
|
See Also |
|
85
|
|
|
-------- |
|
86
|
|
|
moist_lapse : Calculate parcel temperature assuming liquid saturation |
|
87
|
|
|
processes |
|
88
|
|
|
parcel_profile : Calculate complete parcel profile |
|
89
|
|
|
potential_temperature |
|
90
|
|
|
|
|
91
|
|
|
""" |
|
92
|
|
|
return temperature * (pressure / pressure[0])**kappa |
|
93
|
|
|
|
|
94
|
|
|
|
|
95
|
|
|
@exporter.export |
|
96
|
|
|
@check_units('[pressure]', '[temperature]') |
|
97
|
|
|
def moist_lapse(pressure, temperature): |
|
98
|
|
|
r"""Calculate the temperature at a level assuming liquid saturation processes. |
|
99
|
|
|
|
|
100
|
|
|
This function lifts a parcel starting at `temperature`. The starting |
|
101
|
|
|
pressure should be the first item in the `pressure` array. Essentially, |
|
102
|
|
|
this function is calculating moist pseudo-adiabats. |
|
103
|
|
|
|
|
104
|
|
|
Parameters |
|
105
|
|
|
---------- |
|
106
|
|
|
pressure : `pint.Quantity` |
|
107
|
|
|
The atmospheric pressure level(s) of interest |
|
108
|
|
|
temperature : `pint.Quantity` |
|
109
|
|
|
The starting temperature |
|
110
|
|
|
|
|
111
|
|
|
Returns |
|
112
|
|
|
------- |
|
113
|
|
|
`pint.Quantity` |
|
114
|
|
|
The temperature corresponding to the the starting temperature and |
|
115
|
|
|
pressure levels. |
|
116
|
|
|
|
|
117
|
|
|
See Also |
|
118
|
|
|
-------- |
|
119
|
|
|
dry_lapse : Calculate parcel temperature assuming dry adiabatic processes |
|
120
|
|
|
parcel_profile : Calculate complete parcel profile |
|
121
|
|
|
|
|
122
|
|
|
Notes |
|
123
|
|
|
----- |
|
124
|
|
|
This function is implemented by integrating the following differential |
|
125
|
|
|
equation: |
|
126
|
|
|
|
|
127
|
|
|
.. math:: \frac{dT}{dP} = \frac{1}{P} \frac{R_d T + L_v r_s} |
|
128
|
|
|
{C_{pd} + \frac{L_v^2 r_s \epsilon}{R_d T^2}} |
|
129
|
|
|
|
|
130
|
|
|
This equation comes from [Bakhshaii2013]_. |
|
131
|
|
|
|
|
132
|
|
|
""" |
|
133
|
|
|
def dt(t, p): |
|
134
|
|
|
t = units.Quantity(t, temperature.units) |
|
135
|
|
|
p = units.Quantity(p, pressure.units) |
|
136
|
|
|
rs = saturation_mixing_ratio(p, t) |
|
137
|
|
|
frac = ((Rd * t + Lv * rs) / |
|
138
|
|
|
(Cp_d + (Lv * Lv * rs * epsilon / (Rd * t * t)))).to('kelvin') |
|
139
|
|
|
return frac / p |
|
140
|
|
|
return units.Quantity(si.odeint(dt, atleast_1d(temperature).squeeze(), |
|
141
|
|
|
pressure.squeeze()).T.squeeze(), temperature.units) |
|
142
|
|
|
|
|
143
|
|
|
|
|
144
|
|
|
@exporter.export |
|
145
|
|
|
@check_units('[pressure]', '[temperature]', '[temperature]') |
|
146
|
|
|
def lcl(pressure, temperature, dewpt, max_iters=50, eps=1e-5): |
|
147
|
|
|
r"""Calculate the lifted condensation level (LCL) using from the starting point. |
|
148
|
|
|
|
|
149
|
|
|
The starting state for the parcel is defined by `temperature`, `dewpt`, |
|
150
|
|
|
and `pressure`. |
|
151
|
|
|
|
|
152
|
|
|
Parameters |
|
153
|
|
|
---------- |
|
154
|
|
|
pressure : `pint.Quantity` |
|
155
|
|
|
The starting atmospheric pressure |
|
156
|
|
|
temperature : `pint.Quantity` |
|
157
|
|
|
The starting temperature |
|
158
|
|
|
dewpt : `pint.Quantity` |
|
159
|
|
|
The starting dew point |
|
160
|
|
|
|
|
161
|
|
|
Returns |
|
162
|
|
|
------- |
|
163
|
|
|
`(pint.Quantity, pint.Quantity)` |
|
164
|
|
|
The LCL pressure and temperature |
|
165
|
|
|
|
|
166
|
|
|
Other Parameters |
|
167
|
|
|
---------------- |
|
168
|
|
|
max_iters : int, optional |
|
169
|
|
|
The maximum number of iterations to use in calculation, defaults to 50. |
|
170
|
|
|
eps : float, optional |
|
171
|
|
|
The desired relative error in the calculated value, defaults to 1e-5. |
|
172
|
|
|
|
|
173
|
|
|
See Also |
|
174
|
|
|
-------- |
|
175
|
|
|
parcel_profile |
|
176
|
|
|
|
|
177
|
|
|
Notes |
|
178
|
|
|
----- |
|
179
|
|
|
This function is implemented using an iterative approach to solve for the |
|
180
|
|
|
LCL. The basic algorithm is: |
|
181
|
|
|
|
|
182
|
|
|
1. Find the dew point from the LCL pressure and starting mixing ratio |
|
183
|
|
|
2. Find the LCL pressure from the starting temperature and dewpoint |
|
184
|
|
|
3. Iterate until convergence |
|
185
|
|
|
|
|
186
|
|
|
The function is guaranteed to finish by virtue of the `max_iters` counter. |
|
187
|
|
|
|
|
188
|
|
|
""" |
|
189
|
|
|
def _lcl_iter(p, p0, w, t): |
|
190
|
|
|
td = dewpoint(vapor_pressure(units.Quantity(p, pressure.units), w)) |
|
191
|
|
|
return (p0 * (td / t) ** (1. / kappa)).m |
|
192
|
|
|
|
|
193
|
|
|
w = mixing_ratio(saturation_vapor_pressure(dewpt), pressure) |
|
194
|
|
|
fp = so.fixed_point(_lcl_iter, pressure.m, args=(pressure.m, w, temperature), |
|
195
|
|
|
xtol=eps, maxiter=max_iters) |
|
196
|
|
|
lcl_p = units.Quantity(fp, pressure.units) |
|
197
|
|
|
return lcl_p, dewpoint(vapor_pressure(lcl_p, w)) |
|
198
|
|
|
|
|
199
|
|
|
|
|
200
|
|
|
@exporter.export |
|
201
|
|
|
@check_units('[pressure]', '[temperature]', '[temperature]') |
|
202
|
|
View Code Duplication |
def lfc(pressure, temperature, dewpt): |
|
|
|
|
|
|
203
|
|
|
r"""Calculate the level of free convection (LFC). |
|
204
|
|
|
|
|
205
|
|
|
This works by finding the first intersection of the ideal parcel path and |
|
206
|
|
|
the measured parcel temperature. |
|
207
|
|
|
|
|
208
|
|
|
Parameters |
|
209
|
|
|
---------- |
|
210
|
|
|
pressure : `pint.Quantity` |
|
211
|
|
|
The atmospheric pressure |
|
212
|
|
|
temperature : `pint.Quantity` |
|
213
|
|
|
The temperature at the levels given by `pressure` |
|
214
|
|
|
dewpt : `pint.Quantity` |
|
215
|
|
|
The dew point at the levels given by `pressure` |
|
216
|
|
|
|
|
217
|
|
|
Returns |
|
218
|
|
|
------- |
|
219
|
|
|
`pint.Quantity` |
|
220
|
|
|
The LFC pressure and temperature |
|
221
|
|
|
|
|
222
|
|
|
See Also |
|
223
|
|
|
-------- |
|
224
|
|
|
parcel_profile |
|
225
|
|
|
|
|
226
|
|
|
""" |
|
227
|
|
|
ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC') |
|
228
|
|
|
|
|
229
|
|
|
# The parcel profile and data have the same first data point, so we ignore |
|
230
|
|
|
# that point to get the real first intersection for the LFC calculation. |
|
231
|
|
|
x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:], |
|
232
|
|
|
direction='increasing') |
|
233
|
|
|
if len(x) == 0: |
|
234
|
|
|
return np.nan * pressure.units, np.nan * temperature.units |
|
235
|
|
|
else: |
|
236
|
|
|
return x[0], y[0] |
|
237
|
|
|
|
|
238
|
|
|
|
|
239
|
|
|
@exporter.export |
|
240
|
|
View Code Duplication |
@check_units('[pressure]', '[temperature]', '[temperature]') |
|
|
|
|
|
|
241
|
|
|
def el(pressure, temperature, dewpt): |
|
242
|
|
|
r"""Calculate the equilibrium level. |
|
243
|
|
|
|
|
244
|
|
|
This works by finding the last intersection of the ideal parcel path and |
|
245
|
|
|
the measured environmental temperature. If there is one or fewer intersections, there is |
|
246
|
|
|
no equilibrium level. |
|
247
|
|
|
|
|
248
|
|
|
Parameters |
|
249
|
|
|
---------- |
|
250
|
|
|
pressure : `pint.Quantity` |
|
251
|
|
|
The atmospheric pressure |
|
252
|
|
|
temperature : `pint.Quantity` |
|
253
|
|
|
The temperature at the levels given by `pressure` |
|
254
|
|
|
dewpt : `pint.Quantity` |
|
255
|
|
|
The dew point at the levels given by `pressure` |
|
256
|
|
|
|
|
257
|
|
|
Returns |
|
258
|
|
|
------- |
|
259
|
|
|
`pint.Quantity, pint.Quantity` |
|
260
|
|
|
The EL pressure and temperature |
|
261
|
|
|
|
|
262
|
|
|
See Also |
|
263
|
|
|
-------- |
|
264
|
|
|
parcel_profile |
|
265
|
|
|
|
|
266
|
|
|
""" |
|
267
|
|
|
ideal_profile = parcel_profile(pressure, temperature[0], dewpt[0]).to('degC') |
|
268
|
|
|
x, y = find_intersections(pressure[1:], ideal_profile[1:], temperature[1:]) |
|
269
|
|
|
|
|
270
|
|
|
# If there is only one intersection, it's the LFC and we return None. |
|
271
|
|
|
if len(x) <= 1: |
|
272
|
|
|
return np.nan * pressure.units, np.nan * temperature.units |
|
273
|
|
|
else: |
|
274
|
|
|
return x[-1], y[-1] |
|
275
|
|
|
|
|
276
|
|
|
|
|
277
|
|
|
@exporter.export |
|
278
|
|
|
@check_units('[pressure]', '[temperature]', '[temperature]') |
|
279
|
|
|
def parcel_profile(pressure, temperature, dewpt): |
|
280
|
|
|
r"""Calculate the profile a parcel takes through the atmosphere. |
|
281
|
|
|
|
|
282
|
|
|
The parcel starts at `temperature`, and `dewpt`, lifted up |
|
283
|
|
|
dry adiabatically to the LCL, and then moist adiabatically from there. |
|
284
|
|
|
`pressure` specifies the pressure levels for the profile. |
|
285
|
|
|
|
|
286
|
|
|
Parameters |
|
287
|
|
|
---------- |
|
288
|
|
|
pressure : `pint.Quantity` |
|
289
|
|
|
The atmospheric pressure level(s) of interest. The first entry should be the starting |
|
290
|
|
|
point pressure. |
|
291
|
|
|
temperature : `pint.Quantity` |
|
292
|
|
|
The starting temperature |
|
293
|
|
|
dewpt : `pint.Quantity` |
|
294
|
|
|
The starting dew point |
|
295
|
|
|
|
|
296
|
|
|
Returns |
|
297
|
|
|
------- |
|
298
|
|
|
`pint.Quantity` |
|
299
|
|
|
The parcel temperatures at the specified pressure levels. |
|
300
|
|
|
|
|
301
|
|
|
See Also |
|
302
|
|
|
-------- |
|
303
|
|
|
lcl, moist_lapse, dry_lapse |
|
304
|
|
|
|
|
305
|
|
|
""" |
|
306
|
|
|
# Find the LCL |
|
307
|
|
|
l = lcl(pressure[0], temperature, dewpt)[0].to(pressure.units) |
|
308
|
|
|
|
|
309
|
|
|
# Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the |
|
310
|
|
|
# LCL is included in the levels. It's slightly redundant in that case, but simplifies |
|
311
|
|
|
# the logic for removing it later. |
|
312
|
|
|
press_lower = concatenate((pressure[pressure >= l], l)) |
|
313
|
|
|
t1 = dry_lapse(press_lower, temperature) |
|
314
|
|
|
|
|
315
|
|
|
# Find moist pseudo-adiabatic profile starting at the LCL |
|
316
|
|
|
press_upper = concatenate((l, pressure[pressure < l])) |
|
317
|
|
|
t2 = moist_lapse(press_upper, t1[-1]).to(t1.units) |
|
318
|
|
|
|
|
319
|
|
|
# Return LCL *without* the LCL point |
|
320
|
|
|
return concatenate((t1[:-1], t2[1:])) |
|
321
|
|
|
|
|
322
|
|
|
|
|
323
|
|
|
@exporter.export |
|
324
|
|
|
@check_units('[pressure]', '[dimensionless]') |
|
325
|
|
|
def vapor_pressure(pressure, mixing): |
|
326
|
|
|
r"""Calculate water vapor (partial) pressure. |
|
327
|
|
|
|
|
328
|
|
|
Given total `pressure` and water vapor `mixing` ratio, calculates the |
|
329
|
|
|
partial pressure of water vapor. |
|
330
|
|
|
|
|
331
|
|
|
Parameters |
|
332
|
|
|
---------- |
|
333
|
|
|
pressure : `pint.Quantity` |
|
334
|
|
|
total atmospheric pressure |
|
335
|
|
|
mixing : `pint.Quantity` |
|
336
|
|
|
dimensionless mass mixing ratio |
|
337
|
|
|
|
|
338
|
|
|
Returns |
|
339
|
|
|
------- |
|
340
|
|
|
`pint.Quantity` |
|
341
|
|
|
The ambient water vapor (partial) pressure in the same units as |
|
342
|
|
|
`pressure`. |
|
343
|
|
|
|
|
344
|
|
|
Notes |
|
345
|
|
|
----- |
|
346
|
|
|
This function is a straightforward implementation of the equation given in many places, |
|
347
|
|
|
such as [Hobbs1977]_ pg.71: |
|
348
|
|
|
|
|
349
|
|
|
.. math:: e = p \frac{r}{r + \epsilon} |
|
350
|
|
|
|
|
351
|
|
|
See Also |
|
352
|
|
|
-------- |
|
353
|
|
|
saturation_vapor_pressure, dewpoint |
|
354
|
|
|
|
|
355
|
|
|
""" |
|
356
|
|
|
return pressure * mixing / (epsilon + mixing) |
|
357
|
|
|
|
|
358
|
|
|
|
|
359
|
|
|
@exporter.export |
|
360
|
|
|
@check_units('[temperature]') |
|
361
|
|
|
def saturation_vapor_pressure(temperature): |
|
362
|
|
|
r"""Calculate the saturation water vapor (partial) pressure. |
|
363
|
|
|
|
|
364
|
|
|
Parameters |
|
365
|
|
|
---------- |
|
366
|
|
|
temperature : `pint.Quantity` |
|
367
|
|
|
The temperature |
|
368
|
|
|
|
|
369
|
|
|
Returns |
|
370
|
|
|
------- |
|
371
|
|
|
`pint.Quantity` |
|
372
|
|
|
The saturation water vapor (partial) pressure |
|
373
|
|
|
|
|
374
|
|
|
See Also |
|
375
|
|
|
-------- |
|
376
|
|
|
vapor_pressure, dewpoint |
|
377
|
|
|
|
|
378
|
|
|
Notes |
|
379
|
|
|
----- |
|
380
|
|
|
Instead of temperature, dewpoint may be used in order to calculate |
|
381
|
|
|
the actual (ambient) water vapor (partial) pressure. |
|
382
|
|
|
|
|
383
|
|
|
The formula used is that from [Bolton1980]_ for T in degrees Celsius: |
|
384
|
|
|
|
|
385
|
|
|
.. math:: 6.112 e^\frac{17.67T}{T + 243.5} |
|
386
|
|
|
|
|
387
|
|
|
""" |
|
388
|
|
|
# Converted from original in terms of C to use kelvin. Using raw absolute values of C in |
|
389
|
|
|
# a formula plays havoc with units support. |
|
390
|
|
|
return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin) / |
|
391
|
|
|
(temperature - 29.65 * units.kelvin)) |
|
392
|
|
|
|
|
393
|
|
|
|
|
394
|
|
|
@exporter.export |
|
395
|
|
|
@check_units('[temperature]', '[dimensionless]') |
|
396
|
|
|
def dewpoint_rh(temperature, rh): |
|
397
|
|
|
r"""Calculate the ambient dewpoint given air temperature and relative humidity. |
|
398
|
|
|
|
|
399
|
|
|
Parameters |
|
400
|
|
|
---------- |
|
401
|
|
|
temperature : `pint.Quantity` |
|
402
|
|
|
Air temperature |
|
403
|
|
|
rh : `pint.Quantity` |
|
404
|
|
|
Relative humidity expressed as a ratio in the range [0, 1] |
|
405
|
|
|
|
|
406
|
|
|
Returns |
|
407
|
|
|
------- |
|
408
|
|
|
`pint.Quantity` |
|
409
|
|
|
The dew point temperature |
|
410
|
|
|
|
|
411
|
|
|
See Also |
|
412
|
|
|
-------- |
|
413
|
|
|
dewpoint, saturation_vapor_pressure |
|
414
|
|
|
|
|
415
|
|
|
""" |
|
416
|
|
|
return dewpoint(rh * saturation_vapor_pressure(temperature)) |
|
417
|
|
|
|
|
418
|
|
|
|
|
419
|
|
|
@exporter.export |
|
420
|
|
|
@check_units('[pressure]') |
|
421
|
|
|
def dewpoint(e): |
|
422
|
|
|
r"""Calculate the ambient dewpoint given the vapor pressure. |
|
423
|
|
|
|
|
424
|
|
|
Parameters |
|
425
|
|
|
---------- |
|
426
|
|
|
e : `pint.Quantity` |
|
427
|
|
|
Water vapor partial pressure |
|
428
|
|
|
|
|
429
|
|
|
Returns |
|
430
|
|
|
------- |
|
431
|
|
|
`pint.Quantity` |
|
432
|
|
|
Dew point temperature |
|
433
|
|
|
|
|
434
|
|
|
See Also |
|
435
|
|
|
-------- |
|
436
|
|
|
dewpoint_rh, saturation_vapor_pressure, vapor_pressure |
|
437
|
|
|
|
|
438
|
|
|
Notes |
|
439
|
|
|
----- |
|
440
|
|
|
This function inverts the [Bolton1980]_ formula for saturation vapor |
|
441
|
|
|
pressure to instead calculate the temperature. This yield the following |
|
442
|
|
|
formula for dewpoint in degrees Celsius: |
|
443
|
|
|
|
|
444
|
|
|
.. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)} |
|
445
|
|
|
|
|
446
|
|
|
""" |
|
447
|
|
|
val = np.log(e / sat_pressure_0c) |
|
448
|
|
|
return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val) |
|
449
|
|
|
|
|
450
|
|
|
|
|
451
|
|
|
@exporter.export |
|
452
|
|
|
@check_units('[pressure]', '[pressure]', '[dimensionless]') |
|
453
|
|
|
def mixing_ratio(part_press, tot_press, molecular_weight_ratio=epsilon): |
|
454
|
|
|
r"""Calculate the mixing ratio of a gas. |
|
455
|
|
|
|
|
456
|
|
|
This calculates mixing ratio given its partial pressure and the total pressure of |
|
457
|
|
|
the air. There are no required units for the input arrays, other than that |
|
458
|
|
|
they have the same units. |
|
459
|
|
|
|
|
460
|
|
|
Parameters |
|
461
|
|
|
---------- |
|
462
|
|
|
part_press : `pint.Quantity` |
|
463
|
|
|
Partial pressure of the constituent gas |
|
464
|
|
|
tot_press : `pint.Quantity` |
|
465
|
|
|
Total air pressure |
|
466
|
|
|
molecular_weight_ratio : `pint.Quantity` or float, optional |
|
467
|
|
|
The ratio of the molecular weight of the constituent gas to that assumed |
|
468
|
|
|
for air. Defaults to the ratio for water vapor to dry air |
|
469
|
|
|
(:math:`\epsilon\approx0.622`). |
|
470
|
|
|
|
|
471
|
|
|
Returns |
|
472
|
|
|
------- |
|
473
|
|
|
`pint.Quantity` |
|
474
|
|
|
The (mass) mixing ratio, dimensionless (e.g. Kg/Kg or g/g) |
|
475
|
|
|
|
|
476
|
|
|
Notes |
|
477
|
|
|
----- |
|
478
|
|
|
This function is a straightforward implementation of the equation given in many places, |
|
479
|
|
|
such as [Hobbs1977]_ pg.73: |
|
480
|
|
|
|
|
481
|
|
|
.. math:: r = \epsilon \frac{e}{p - e} |
|
482
|
|
|
|
|
483
|
|
|
See Also |
|
484
|
|
|
-------- |
|
485
|
|
|
saturation_mixing_ratio, vapor_pressure |
|
486
|
|
|
|
|
487
|
|
|
""" |
|
488
|
|
|
return molecular_weight_ratio * part_press / (tot_press - part_press) |
|
489
|
|
|
|
|
490
|
|
|
|
|
491
|
|
|
@exporter.export |
|
492
|
|
|
@check_units('[pressure]', '[temperature]') |
|
493
|
|
|
def saturation_mixing_ratio(tot_press, temperature): |
|
494
|
|
|
r"""Calculate the saturation mixing ratio of water vapor. |
|
495
|
|
|
|
|
496
|
|
|
This calculation is given total pressure and the temperature. The implementation |
|
497
|
|
|
uses the formula outlined in [Hobbs1977]_ pg.73. |
|
498
|
|
|
|
|
499
|
|
|
Parameters |
|
500
|
|
|
---------- |
|
501
|
|
|
tot_press: `pint.Quantity` |
|
502
|
|
|
Total atmospheric pressure |
|
503
|
|
|
temperature: `pint.Quantity` |
|
504
|
|
|
The temperature |
|
505
|
|
|
|
|
506
|
|
|
Returns |
|
507
|
|
|
------- |
|
508
|
|
|
`pint.Quantity` |
|
509
|
|
|
The saturation mixing ratio, dimensionless |
|
510
|
|
|
|
|
511
|
|
|
""" |
|
512
|
|
|
return mixing_ratio(saturation_vapor_pressure(temperature), tot_press) |
|
513
|
|
|
|
|
514
|
|
|
|
|
515
|
|
|
@exporter.export |
|
516
|
|
|
@check_units('[pressure]', '[temperature]') |
|
517
|
|
|
def equivalent_potential_temperature(pressure, temperature): |
|
518
|
|
|
r"""Calculate equivalent potential temperature. |
|
519
|
|
|
|
|
520
|
|
|
This calculation must be given an air parcel's pressure and temperature. |
|
521
|
|
|
The implementation uses the formula outlined in [Hobbs1977]_ pg.78-79. |
|
522
|
|
|
|
|
523
|
|
|
Parameters |
|
524
|
|
|
---------- |
|
525
|
|
|
pressure: `pint.Quantity` |
|
526
|
|
|
Total atmospheric pressure |
|
527
|
|
|
temperature: `pint.Quantity` |
|
528
|
|
|
The temperature |
|
529
|
|
|
|
|
530
|
|
|
Returns |
|
531
|
|
|
------- |
|
532
|
|
|
`pint.Quantity` |
|
533
|
|
|
The corresponding equivalent potential temperature of the parcel |
|
534
|
|
|
|
|
535
|
|
|
Notes |
|
536
|
|
|
----- |
|
537
|
|
|
.. math:: \Theta_e = \Theta e^\frac{L_v r_s}{C_{pd} T} |
|
538
|
|
|
|
|
539
|
|
|
""" |
|
540
|
|
|
pottemp = potential_temperature(pressure, temperature) |
|
541
|
|
|
smixr = saturation_mixing_ratio(pressure, temperature) |
|
542
|
|
|
return pottemp * np.exp(Lv * smixr / (Cp_d * temperature)) |
|
543
|
|
|
|
|
544
|
|
|
|
|
545
|
|
|
@exporter.export |
|
546
|
|
|
@check_units('[temperature]', '[dimensionless]', '[dimensionless]') |
|
547
|
|
|
def virtual_temperature(temperature, mixing, molecular_weight_ratio=epsilon): |
|
548
|
|
|
r"""Calculate virtual temperature. |
|
549
|
|
|
|
|
550
|
|
|
This calculation must be given an air parcel's temperature and mixing ratio. |
|
551
|
|
|
The implementation uses the formula outlined in [Hobbs2006]_ pg.80. |
|
552
|
|
|
|
|
553
|
|
|
Parameters |
|
554
|
|
|
---------- |
|
555
|
|
|
temperature: `pint.Quantity` |
|
556
|
|
|
The temperature |
|
557
|
|
|
mixing : `pint.Quantity` |
|
558
|
|
|
dimensionless mass mixing ratio |
|
559
|
|
|
molecular_weight_ratio : `pint.Quantity` or float, optional |
|
560
|
|
|
The ratio of the molecular weight of the constituent gas to that assumed |
|
561
|
|
|
for air. Defaults to the ratio for water vapor to dry air |
|
562
|
|
|
(:math:`\epsilon\approx0.622`). |
|
563
|
|
|
|
|
564
|
|
|
Returns |
|
565
|
|
|
------- |
|
566
|
|
|
`pint.Quantity` |
|
567
|
|
|
The corresponding virtual temperature of the parcel |
|
568
|
|
|
|
|
569
|
|
|
Notes |
|
570
|
|
|
----- |
|
571
|
|
|
.. math:: T_v = T \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})} |
|
572
|
|
|
|
|
573
|
|
|
""" |
|
574
|
|
|
return temperature * ((mixing + molecular_weight_ratio) / |
|
575
|
|
|
(molecular_weight_ratio * (1 + mixing))) |
|
576
|
|
|
|
|
577
|
|
|
|
|
578
|
|
|
@exporter.export |
|
579
|
|
|
@check_units('[pressure]', '[temperature]', '[dimensionless]', '[dimensionless]') |
|
580
|
|
|
def virtual_potential_temperature(pressure, temperature, mixing, |
|
581
|
|
|
molecular_weight_ratio=epsilon): |
|
582
|
|
|
r"""Calculate virtual potential temperature. |
|
583
|
|
|
|
|
584
|
|
|
This calculation must be given an air parcel's pressure, temperature, and mixing ratio. |
|
585
|
|
|
The implementation uses the formula outlined in [Markowski2010]_ pg.13. |
|
586
|
|
|
|
|
587
|
|
|
Parameters |
|
588
|
|
|
---------- |
|
589
|
|
|
pressure: `pint.Quantity` |
|
590
|
|
|
Total atmospheric pressure |
|
591
|
|
|
temperature: `pint.Quantity` |
|
592
|
|
|
The temperature |
|
593
|
|
|
mixing : `pint.Quantity` |
|
594
|
|
|
dimensionless mass mixing ratio |
|
595
|
|
|
molecular_weight_ratio : `pint.Quantity` or float, optional |
|
596
|
|
|
The ratio of the molecular weight of the constituent gas to that assumed |
|
597
|
|
|
for air. Defaults to the ratio for water vapor to dry air |
|
598
|
|
|
(:math:`\epsilon\approx0.622`). |
|
599
|
|
|
|
|
600
|
|
|
Returns |
|
601
|
|
|
------- |
|
602
|
|
|
`pint.Quantity` |
|
603
|
|
|
The corresponding virtual potential temperature of the parcel |
|
604
|
|
|
|
|
605
|
|
|
Notes |
|
606
|
|
|
----- |
|
607
|
|
|
.. math:: \Theta_v = \Theta \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})} |
|
608
|
|
|
|
|
609
|
|
|
""" |
|
610
|
|
|
pottemp = potential_temperature(pressure, temperature) |
|
611
|
|
|
return virtual_temperature(pottemp, mixing, molecular_weight_ratio) |
|
612
|
|
|
|
|
613
|
|
|
|
|
614
|
|
|
@exporter.export |
|
615
|
|
|
@check_units('[pressure]', '[temperature]', '[dimensionless]', '[dimensionless]') |
|
616
|
|
|
def density(pressure, temperature, mixing, molecular_weight_ratio=epsilon): |
|
617
|
|
|
r"""Calculate density. |
|
618
|
|
|
|
|
619
|
|
|
This calculation must be given an air parcel's pressure, temperature, and mixing ratio. |
|
620
|
|
|
The implementation uses the formula outlined in [Hobbs2006]_ pg.67. |
|
621
|
|
|
|
|
622
|
|
|
Parameters |
|
623
|
|
|
---------- |
|
624
|
|
|
temperature: `pint.Quantity` |
|
625
|
|
|
The temperature |
|
626
|
|
|
pressure: `pint.Quantity` |
|
627
|
|
|
Total atmospheric pressure |
|
628
|
|
|
mixing : `pint.Quantity` |
|
629
|
|
|
dimensionless mass mixing ratio |
|
630
|
|
|
molecular_weight_ratio : `pint.Quantity` or float, optional |
|
631
|
|
|
The ratio of the molecular weight of the constituent gas to that assumed |
|
632
|
|
|
for air. Defaults to the ratio for water vapor to dry air |
|
633
|
|
|
(:math:`\epsilon\approx0.622`). |
|
634
|
|
|
|
|
635
|
|
|
Returns |
|
636
|
|
|
------- |
|
637
|
|
|
`pint.Quantity` |
|
638
|
|
|
The corresponding density of the parcel |
|
639
|
|
|
|
|
640
|
|
|
Notes |
|
641
|
|
|
----- |
|
642
|
|
|
.. math:: \rho = \frac{p}{R_dT_v} |
|
643
|
|
|
|
|
644
|
|
|
""" |
|
645
|
|
|
virttemp = virtual_temperature(temperature, mixing, molecular_weight_ratio) |
|
646
|
|
|
return (pressure / (Rd * virttemp)).to(units.kilogram / units.meter ** 3) |
|
647
|
|
|
|
|
648
|
|
|
|
|
649
|
|
|
@exporter.export |
|
650
|
|
|
@check_units('[temperature]', '[temperature]', '[pressure]') |
|
651
|
|
|
def relative_humidity_wet_psychrometric(dry_bulb_temperature, web_bulb_temperature, |
|
652
|
|
|
pressure, **kwargs): |
|
653
|
|
|
r"""Calculate the relative humidity with wet bulb and dry bulb temperatures. |
|
654
|
|
|
|
|
655
|
|
|
This uses a psychrometric relationship as outlined in [WMO8-2008]_, with |
|
656
|
|
|
coefficients from [Fan1987]_. |
|
657
|
|
|
|
|
658
|
|
|
Parameters |
|
659
|
|
|
---------- |
|
660
|
|
|
dry_bulb_temperature: `pint.Quantity` |
|
661
|
|
|
Dry bulb temperature |
|
662
|
|
|
web_bulb_temperature: `pint.Quantity` |
|
663
|
|
|
Wet bulb temperature |
|
664
|
|
|
pressure: `pint.Quantity` |
|
665
|
|
|
Total atmospheric pressure |
|
666
|
|
|
|
|
667
|
|
|
Returns |
|
668
|
|
|
------- |
|
669
|
|
|
`pint.Quantity` |
|
670
|
|
|
Relative humidity |
|
671
|
|
|
|
|
672
|
|
|
Notes |
|
673
|
|
|
----- |
|
674
|
|
|
.. math:: RH = 100 \frac{e}{e_s} |
|
675
|
|
|
|
|
676
|
|
|
* :math:`RH` is relative humidity |
|
677
|
|
|
* :math:`e` is vapor pressure from the wet psychrometric calculation |
|
678
|
|
|
* :math:`e_s` is the saturation vapor pressure |
|
679
|
|
|
|
|
680
|
|
|
See Also |
|
681
|
|
|
-------- |
|
682
|
|
|
psychrometric_vapor_pressure_wet, saturation_vapor_pressure |
|
683
|
|
|
|
|
684
|
|
|
""" |
|
685
|
|
|
return (100 * units.percent * psychrometric_vapor_pressure_wet(dry_bulb_temperature, |
|
686
|
|
|
web_bulb_temperature, pressure, **kwargs) / |
|
687
|
|
|
saturation_vapor_pressure(dry_bulb_temperature)) |
|
688
|
|
|
|
|
689
|
|
|
|
|
690
|
|
|
@exporter.export |
|
691
|
|
|
@check_units('[temperature]', '[temperature]', '[pressure]') |
|
692
|
|
|
def psychrometric_vapor_pressure_wet(dry_bulb_temperature, wet_bulb_temperature, pressure, |
|
693
|
|
|
psychrometer_coefficient=6.21e-4 / units.kelvin): |
|
694
|
|
|
r"""Calculate the vapor pressure with wet bulb and dry bulb temperatures. |
|
695
|
|
|
|
|
696
|
|
|
This uses a psychrometric relationship as outlined in [WMO8-2008]_, with |
|
697
|
|
|
coefficients from [Fan1987]_. |
|
698
|
|
|
|
|
699
|
|
|
Parameters |
|
700
|
|
|
---------- |
|
701
|
|
|
dry_bulb_temperature: `pint.Quantity` |
|
702
|
|
|
Dry bulb temperature |
|
703
|
|
|
wet_bulb_temperature: `pint.Quantity` |
|
704
|
|
|
Wet bulb temperature |
|
705
|
|
|
pressure: `pint.Quantity` |
|
706
|
|
|
Total atmospheric pressure |
|
707
|
|
|
psychrometer_coefficient: `pint.Quantity` |
|
708
|
|
|
Psychrometer coefficient |
|
709
|
|
|
|
|
710
|
|
|
Returns |
|
711
|
|
|
------- |
|
712
|
|
|
`pint.Quantity` |
|
713
|
|
|
Vapor pressure |
|
714
|
|
|
|
|
715
|
|
|
Notes |
|
716
|
|
|
----- |
|
717
|
|
|
.. math:: e' = e'_w(T_w) - A p (T - T_w) |
|
718
|
|
|
|
|
719
|
|
|
* :math:`e'` is vapor pressure |
|
720
|
|
|
* :math:`e'_w(T_w)` is the saturation vapor pressure with respect to water at temperature |
|
721
|
|
|
:math:`T_w` |
|
722
|
|
|
* :math:`p` is the pressure of the wet bulb |
|
723
|
|
|
* :math:`T` is the temperature of the dry bulb |
|
724
|
|
|
* :math:`T_w` is the temperature of the wet bulb |
|
725
|
|
|
* :math:`A` is the psychrometer coefficient |
|
726
|
|
|
|
|
727
|
|
|
Psychrometer coefficient depends on the specific instrument being used and the ventilation |
|
728
|
|
|
of the instrument. |
|
729
|
|
|
|
|
730
|
|
|
See Also |
|
731
|
|
|
-------- |
|
732
|
|
|
saturation_vapor_pressure |
|
733
|
|
|
|
|
734
|
|
|
""" |
|
735
|
|
|
return (saturation_vapor_pressure(wet_bulb_temperature) - psychrometer_coefficient * |
|
736
|
|
|
pressure * (dry_bulb_temperature - wet_bulb_temperature).to('kelvin')) |
|
737
|
|
|
|
|
738
|
|
|
|
|
739
|
|
|
@exporter.export |
|
740
|
|
|
@check_units('[dimensionless]', '[temperature]', '[pressure]') |
|
741
|
|
|
def relative_humidity_from_mixing_ratio(mixing_ratio, temperature, pressure): |
|
742
|
|
|
r"""Calculate the relative humidity from mixing ratio, temperature, and pressure. |
|
743
|
|
|
|
|
744
|
|
|
Parameters |
|
745
|
|
|
---------- |
|
746
|
|
|
mixing_ratio: `pint.Quantity` |
|
747
|
|
|
Dimensionless mass mixing ratio |
|
748
|
|
|
temperature: `pint.Quantity` |
|
749
|
|
|
Air temperature |
|
750
|
|
|
pressure: `pint.Quantity` |
|
751
|
|
|
Total atmospheric pressure |
|
752
|
|
|
|
|
753
|
|
|
Returns |
|
754
|
|
|
------- |
|
755
|
|
|
`pint.Quantity` |
|
756
|
|
|
Relative humidity |
|
757
|
|
|
|
|
758
|
|
|
Notes |
|
759
|
|
|
----- |
|
760
|
|
|
Formula from [Hobbs 1977]_ pg. 74. |
|
761
|
|
|
.. math:: RH = 100 \frac{w}{w_s} |
|
762
|
|
|
|
|
763
|
|
|
* :math:`RH` is relative humidity |
|
764
|
|
|
* :math:`w` is mxing ratio |
|
765
|
|
|
* :math:`w_s` is the saturation mixing ratio |
|
766
|
|
|
|
|
767
|
|
|
See Also |
|
768
|
|
|
-------- |
|
769
|
|
|
saturation_mixing_ratio |
|
770
|
|
|
|
|
771
|
|
|
""" |
|
772
|
|
|
return (100 * units.percent * |
|
773
|
|
|
mixing_ratio / saturation_mixing_ratio(pressure, temperature)) |
|
774
|
|
|
|
|
775
|
|
|
|
|
776
|
|
|
@exporter.export |
|
777
|
|
|
@check_units('[dimensionless]') |
|
778
|
|
|
def mixing_ratio_from_specific_humidity(specific_humidity): |
|
779
|
|
|
r"""Calculate the mixing ratio from specific humidity. |
|
780
|
|
|
|
|
781
|
|
|
Parameters |
|
782
|
|
|
---------- |
|
783
|
|
|
specific_humidity: `pint.Quantity` |
|
784
|
|
|
Specific humidity of air |
|
785
|
|
|
|
|
786
|
|
|
Returns |
|
787
|
|
|
------- |
|
788
|
|
|
`pint.Quantity` |
|
789
|
|
|
Mixing ratio |
|
790
|
|
|
|
|
791
|
|
|
Notes |
|
792
|
|
|
----- |
|
793
|
|
|
Formula from [Salby1996]_ pg. 118. |
|
794
|
|
|
.. math:: w = \frac{q}{1-q} |
|
795
|
|
|
|
|
796
|
|
|
* :math:`w` is mxing ratio |
|
797
|
|
|
* :math:`q` is the specific humidity |
|
798
|
|
|
|
|
799
|
|
|
See Also |
|
800
|
|
|
-------- |
|
801
|
|
|
mixing_ratio |
|
802
|
|
|
|
|
803
|
|
|
""" |
|
804
|
|
|
return specific_humidity / (1 - specific_humidity) |
|
805
|
|
|
|
|
806
|
|
|
|
|
807
|
|
|
@exporter.export |
|
808
|
|
|
@check_units('[dimensionless]', '[temperature]', '[pressure]') |
|
809
|
|
|
def relative_humidity_from_specific_humidity(specific_humidity, temperature, pressure): |
|
810
|
|
|
r"""Calculate the relative humidity from specific humidity, temperature, and pressure. |
|
811
|
|
|
|
|
812
|
|
|
Parameters |
|
813
|
|
|
---------- |
|
814
|
|
|
specific_humidity: `pint.Quantity` |
|
815
|
|
|
Specific humidity of air |
|
816
|
|
|
temperature: `pint.Quantity` |
|
817
|
|
|
Air temperature |
|
818
|
|
|
pressure: `pint.Quantity` |
|
819
|
|
|
Total atmospheric pressure |
|
820
|
|
|
|
|
821
|
|
|
Returns |
|
822
|
|
|
------- |
|
823
|
|
|
`pint.Quantity` |
|
824
|
|
|
Relative humidity |
|
825
|
|
|
|
|
826
|
|
|
Notes |
|
827
|
|
|
----- |
|
828
|
|
|
Formula from [Hobbs 1977]_ pg. 74. and [Salby1996]_ pg. 118. |
|
829
|
|
|
.. math:: RH = 100 \frac{q}{(1-q)w_s} |
|
830
|
|
|
|
|
831
|
|
|
* :math:`RH` is relative humidity |
|
832
|
|
|
* :math:`q` is specific humidity |
|
833
|
|
|
* :math:`w_s` is the saturation mixing ratio |
|
834
|
|
|
|
|
835
|
|
|
See Also |
|
836
|
|
|
-------- |
|
837
|
|
|
relative_humidity_from_mixing_ratio |
|
838
|
|
|
|
|
839
|
|
|
""" |
|
840
|
|
|
return (100 * units.percent * |
|
841
|
|
|
mixing_ratio_from_specific_humidity(specific_humidity) / |
|
842
|
|
|
saturation_mixing_ratio(pressure, temperature)) |
|
843
|
|
|
|