import numpy as np
from scipy.optimize import brentq
from scipy.interpolate import interp1d
from derive_functions import qv_from_vapr, temp_from_vapr
import iterative
[docs]def compute_lcl(sfc_T, sfc_p, sfc_r=None, sfc_td=None):
"""
Compute the LCL pressure and temperature given parcel temperature, pressure, and mixing ratio or dewpoint.
Args:
sfc_T: Parcel temperature in degrees C
sfc_p: Parcel pressure in hPa
sfc_r: Parcel mixing ratio in g/g. Optional if sfc_td is specified.
sfc_td: Parcel dewpoint in degrees C. Optional if sfc_r is specified.
Returns:
LCL temperature in degrees C and LCL pressure in hPa.
"""
if sfc_r is None and sfc_td is None:
print "computeLCL() must receive either sfc_r or sfc_td"
return np.nan, np.nan
if sfc_r is None:
sfc_r = qv_from_vapr(t=(sfc_td + 273.15), p=(sfc_p * 100))
lcl_T, lcl_p, lcl_r = iterative.solve_iteratively(iterative.condensation_temp, [sfc_T + 273.15, sfc_p / 10., sfc_r])
return lcl_T - 273.15, lcl_p * 10
[docs]def compute_lfc(parcel, snd):
"""
Computes LFC given a parcel trace and environmental trace.
Args:
parcel (tuple): A tuple of (parcel pressure in hPa, parcel temperature in degrees C). Both should be 1-
dimensional arrays.
snd (tuple): A tuple of (environmental pressure in hPa, environmental temperature in degrees C). Both should be
1-dimensional arrays.
Returns:
The pressure of the parcel LFC in hPa. In the case of multiple LFCs, the highest one (lowest in pressure) is
returned.
"""
pres, (t_parcel, t_snd) = _merge_profiles(parcel, snd)
equil = _find_equilibria(pres, (t_parcel, t_snd), direction='np')
if len(equil) > 0:
lfc = equil[-1]
else:
lfc = np.nan
return lfc
[docs]def compute_el(parcel, snd):
"""
Computes EL given a parcel trace and environmental trace.
Args:
parcel (tuple): A tuple of (parcel pressure in hPa, parcel temperature in degrees C). Both should be 1-
dimensional arrays.
snd (tuple): A tuple of (environmental pressure in hPa, environmental temperature in degrees C). Both should be
1-dimensional arrays.
Returns:
The pressure of the parcel EL in hPa. In the case of multiple ELs, the highest one (lowest in pressure) is
returned.
"""
pres, (t_parcel, t_snd) = _merge_profiles(parcel, snd)
equil = _find_equilibria(pres, (t_parcel, t_snd), direction='pn')
if len(equil) > 0:
el = equil[-1]
else:
el = np.nan
return el
[docs]def sb_parcel(t_snd, p_snd, td_snd):
"""
Find the surface-based parcel given temperature, dewpoint, and pressure profiles.
Args:
t_snd (np.ndarray): A 1-dimensional array of temperatures in degrees C.
p_snd (np.ndarray): A 1-dimensional array of pressures in hPa.
td_snd (np.ndarray): A 1-dimensional array of dewpoints in degrees C.
Returns:
Temperature in degrees C, pressure in hPa, and dewpoint in degrees C of the surface-based parcel.
"""
return t_snd[0], p_snd[0], td_snd[0]
[docs]def ml_parcel(t_snd, p_snd, td_snd, ml_depth=100):
"""
Find the mixed-layer parcel given temperature, dewpoint, and pressure profiles.
Args:
t_snd (np.ndarray): A 1-dimensional array of temperatures in degrees C.
p_snd (np.ndarray): A 1-dimensional array of pressures in hPa.
td_snd (np.ndarray): A 1-dimensional array of dewpoints in degrees C.
ml_depth (float): The depth of the mixed layer in hPa. Optional, defaults to 100 hPa.
Returns:
Temperature, in degrees C, pressure in hPa, and dewpoint in degrees C of the mixed-layer parcel.
"""
if ml_depth <= 0:
raise ValueError("ml_parcel must have a mixed layer depth > 0")
r_snd = qv_from_vapr(t=(t_snd + 273.15), p=(p_snd * 100))
pt_snd = (t_snd + 273.15) * (1000. / p_snd) ** (2. / 7.)
ml_idx = np.where(p_snd >= (p_snd[0] - ml_depth))
ml_pt = (pt_snd[ml_idx] * p_snd[ml_idx]).sum() / p_snd[ml_idx].sum()
ml_r = (r_snd[ml_idx] * p_snd[ml_idx]).sum() / p_snd[ml_idx].sum()
t_prc = ml_pt * (p_snd[0] / 1000.) ** (2. / 7.) - 273.15
td_prc = temp_from_vapr(qv=ml_r, p=(p_snd[0] * 100.))
return t_prc, p_snd[0], td_prc
[docs]def mu_parcel(t_snd, p_snd, td_snd, search_depth=300):
"""
Find the most-unstable parcel given temperature, dewpoint, and pressure profiles.
Args:
t_snd (np.ndarray): A 1-dimensional array of temperatures in degrees C.
p_snd (np.ndarray): A 1-dimensional array of pressures in hPa.
td_snd (np.ndarray): A 1-dimensional array of dewpoints in degrees C.
search_depth (float): The depth of the layer to search in hPa. Optional, defaults to 300 hPa.
Returns:
Temperature, in degrees C, pressure in hPa, and dewpoint in degrees C of the most-unstable parcel.
"""
return
def _merge_profiles(*profs):
pres_profs, quant_profs = zip(*profs)
pres_merge = np.unique(np.concatenate(pres_profs))[::-1]
quant_merge = []
for pprof, qprof in profs:
qmerge = interp1d(pprof[::-1], qprof[::-1], bounds_error=False, fill_value=np.nan)(pres_merge)
quant_merge.append(qmerge)
return pres_merge, quant_merge
def _find_equilibria(pres, temp, direction='np'):
t_parcel, t_snd = temp
t_diff = t_parcel - t_snd
t_crossings = np.where((t_diff[:-1] * t_diff[1:]) < 0)[0]
if direction == 'np':
t_cross_keep = np.array([ idx for idx in t_crossings if t_diff[idx] < 0 ])
elif direction == 'pn':
t_cross_keep = np.array([ idx for idx in t_crossings if t_diff[idx] > 0 ])
else:
raise ValueError("_find_equilibria() got '%s' for direction. Should be one of 'np' and 'pn'." % direction)
tdiff_interp = interp1d(pres[::-1], t_diff[::-1])
equil = []
for tck in t_cross_keep:
eql = brentq(tdiff_interp, pres[tck], pres[tck + 1])
equil.append(eql)
return equil