diff --git a/lib/xmet/thermo.py b/lib/xmet/thermo.py index e96d02a..d1e3758 100644 --- a/lib/xmet/thermo.py +++ b/lib/xmet/thermo.py @@ -1,6 +1,12 @@ +from typing import Callable + LAPSE_RATE_DRY = 9.8 / 1000 # degrees C per km LAPSE_RATE_MOIST = 4.0 / 1000 +PRESSURE_MAX = 1050 # millibar +PRESSURE_MIN = 100 +PRESSURE_STEP = 50 + def kelvin(c: float) -> float: return 273.15 + c @@ -67,3 +73,45 @@ def moist_lapse_rate(temp: float, dewpoint: float, pressure: float) -> float: return g * (1 + (Hv * r) / (Rsd * T)) \ / (Cpd + (((Hv**2) * r) / (Rsw * (T**2)))) + +def loft_parcel(start_temp: float, + start_pressure: float, + lapse_rate: Callable): + """ + Loft a parcel of air from a given pressure, at a given temperature, + yielding a Tuple containing the temperature and pressure of that parcel + at each increment of elevation. A Callable which, given a temperature, + dewpoint, pressure value, will be dispatched to obtain the lapse rate + for each height. + """ + start_height = pressure_height(start_pressure) + sx_last = None + sy_last = None + height_last = None + + temp = start_temp + pressure = start_pressure + + while pressure >= PRESSURE_MIN: + height = pressure_height(pressure) + + if height_last is None: + height_last = height + + try: + rate = lapse_rate(temp, temp, pressure) + except OverflowError: + break + + temp = lapse(temp, height - height_last, rate) + + yield temp, pressure + + height_last = height + + if pressure == PRESSURE_MIN: + break + elif pressure - 10.0 < PRESSURE_MIN: + pressure = PRESSURE_MIN + else: + pressure -= 10.0