Adds additional lapse rate options for moist_lapse

src/metpy/calc/thermo.py

@@ -266,7 +266,7 @@ def dry_lapse(pressure, temperature, reference_pressure=None, vertical_dim=0):
     },
     '[temperature]'
 )
-def moist_lapse(pressure, temperature, reference_pressure=None):
+def moist_lapse(pressure, temperature, reference_pressure=None, lapse_type='standard', params=None):
     r"""Calculate the temperature at a level assuming liquid saturation processes.
 
     This function lifts a parcel starting at `temperature`. The starting pressure can
@@ -285,6 +285,29 @@ def moist_lapse(pressure, temperature, reference_pressure=None):
         Reference pressure; if not given, it defaults to the first element of the
         pressure array.
 
+    lapse_type : `string`, optional
+        Definition of moist adiabat to use; if not given, it defaults to moist_lapse
+        Options:
+            'standard' for simplified pseudoadiabatic process
+            'pseudoadiabatic' for pseudoadiabatic moist process
+            'reversible' for reversible moist process
+            'so13' for Singh and O'Gorman (2013);  https://doi.org/10.1002/grl.50796
+            'r14' for Romps (2014); https://doi.org/10.1175/JCLI-D-14-00255.1
+        More info: https://glossary.ametsoc.org/wiki/Adiabatic_lapse_rate
+
+    params : `dict` or None, optional
+        External parameters used for the some lapse_types
+        Required parameters:
+            For 'so13': {
+                'ep0': scalar, entrainment constant [unitless],
+                'rh0': scalar, ambient relative humidity [unitless],
+                }
+            For 'r14': {
+                'de': scalar or 1-d array, detrainment rate [m**-1],
+                'ep': scalar or 1-d array, entrainment rate [m**-1],
+                'pa': 1-d array, optional, pressure levels defining detrainment and entrainment profile [Pa]
+                }
+
     Returns
     -------
     `pint.Quantity`
@@ -321,22 +344,79 @@ def moist_lapse(pressure, temperature, reference_pressure=None):
        Renamed ``ref_pressure`` parameter to ``reference_pressure``
 
     """
-    def dt(p, t):
+    def dt_standard(p, t, params):
         rs = saturation_mixing_ratio._nounit(p, t)
         frac = (
             (mpconsts.nounit.Rd * t + mpconsts.nounit.Lv * rs)
             / (mpconsts.nounit.Cp_d + (
-                mpconsts.nounit.Lv * mpconsts.nounit.Lv * rs * mpconsts.nounit.epsilon
+                mpconsts.nounit.Lv**2 * rs * mpconsts.nounit.epsilon
+                / (mpconsts.nounit.Rd * t**2)
+            ))
+        )
+        return frac / p
+
+    def dt_pseudoadiabatic(p, t, params):
+        rs = saturation_mixing_ratio._nounit(p, t)
+        frac = ( (1 + rs)*(mpconsts.nounit.Rd * t + mpconsts.nounit.Lv * rs)
+                / (mpconsts.nounit.Cp_d + rs*mpconsts.nounit.Cv_d +  (mpconsts.nounit.Lv**2 * rs * (mpconsts.nounit.epsilon + rs)
+                                    / (mpconsts.nounit.Rd * t**2))))
+        return frac / p
+
+    def dt_reversible(p, t, params):
+        rs = saturation_mixing_ratio._nounit(p, t)
+        rl = params['rt'] - rs  ## assuming no ice content
+        frac = ( (1 + params['rt'])*(mpconsts.nounit.Rd * t + mpconsts.nounit.Lv * rs)
+                / (mpconsts.nounit.Cp_d + rs*mpconsts.nounit.Cv_d + rl*mpconsts.nounit.Cp_l + (mpconsts.nounit.Lv**2 * rs * (mpconsts.nounit.epsilon + rs)
+                                    / (mpconsts.nounit.Rd * t**2))))
+        return frac / p
+
+    def dt_so13(p, t, params):
+        zp = -params['h0']*np.log(p/params['p0']) # pseudoheight
+        ep = params['ep0']/zp # entrainment rate
+        rs = saturation_mixing_ratio._nounit(p, t)
+        qs = specific_humidity_from_mixing_ratio(rs)
+        frac = (
+            (mpconsts.nounit.Rd*t + mpconsts.nounit.Lv*qs + ep*qs*mpconsts.nounit.Lv*(1-params['rh0'])*mpconsts.nounit.Rd*t/mpconsts.nounit.g)
+            / (mpconsts.nounit.Cp_d + (
+                mpconsts.nounit.Lv**2 * qs * mpconsts.nounit.epsilon
                 / (mpconsts.nounit.Rd * t**2)
             ))
         )
         return frac / p
 
+    def dt_r14(p, t, params):
+        ep = np.interp(p,params['pa'],params['ep']) if hasattr(params['ep'],'__len__') else params['ep'] # entrainment rate at p
+        de = np.interp(p,params['pa'],params['de']) if hasattr(params['de'],'__len__') else params['de'] # detrainment rate at p
+        rs = saturation_mixing_ratio._nounit(p, t)
+        qs = specific_humidity_from_mixing_ratio(rs)
+        a1 = mpconsts.nounit.Rv*mpconsts.nounit.Cp_d*t**2/mpconsts.nounit.Lv + qs*mpconsts.nounit.Lv
+        a2 = mpconsts.nounit.Rv*mpconsts.nounit.Cp_d*t**2/mpconsts.nounit.Lv*(de+mpconsts.nounit.g/(mpconsts.nounit.Rd*t)) + qs*mpconsts.nounit.Lv*(de-ep) - mpconsts.nounit.g
+        a3 = (mpconsts.nounit.Rv*mpconsts.nounit.Cp_d*t/(mpconsts.nounit.Rd*mpconsts.nounit.Lv) - 1)*mpconsts.nounit.g*de
+        frac = mpconsts.nounit.Rd*t/(mpconsts.nounit.g) * mpconsts.nounit.Rv*t**2/mpconsts.nounit.Lv * ((-a2+np.sqrt(a2**2-4*a1*a3))/(2*a1) + mpconsts.nounit.g/(mpconsts.nounit.Rd*t))
+        return frac / p
+
     temperature = np.atleast_1d(temperature)
     pressure = np.atleast_1d(pressure)
     if reference_pressure is None:
         reference_pressure = pressure[0]
 
+    if lapse_type == 'standard':
+        dt=dt_standard
+    elif lapse_type == 'pseudoadiabatic':
+        dt=dt_pseudoadiabatic
+    elif lapse_type == 'reversible':
+        dt=dt_reversible
+        params={'rt':saturation_mixing_ratio._nounit(reference_pressure,temperature)} # total water at LCL = rs
+    elif lapse_type == 'so13':
+        dt=dt_so13
+        params.update{{'h0':mpconsts.nounit.Rd*temperature[0]/mpconsts.nounit.g, 'p0':pressure[0]}}
+    elif lapse_type == 'r14':
+        dt=dt_r14
+    else:
+        raise ValueError('Specified lapse_type is not supported. '
+                         'Choose from standard, pseudoadiabatic, reversible, '
+                         'so13, or r14.')
+
     if np.isnan(reference_pressure) or np.all(np.isnan(temperature)):
         return np.full((temperature.size, pressure.size), np.nan)
 
@@ -347,7 +427,7 @@ def dt(p, t):
 
     # It would be preferable to use a regular solver like RK45, but as of scipy 1.8.0
     # anything other than LSODA goes into an infinite loop when given NaNs for y0.
-    solver_args = {'fun': dt, 'y0': temperature,
+    solver_args = {'fun':lambda p,t:dt(p,t,params), 'y0': temperature,
                    'method': 'LSODA', 'atol': 1e-7, 'rtol': 1.5e-8}
 
     # Need to handle close points to avoid an error in the solver
@@ -911,11 +991,10 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
         return (units.Quantity(np.nan, pressure.units),
                 units.Quantity(np.nan, temperature.units))
 
-
 @exporter.export
 @preprocess_and_wrap(wrap_like='pressure')
 @check_units('[pressure]', '[temperature]', '[temperature]')
-def parcel_profile(pressure, temperature, dewpoint):
+def parcel_profile(pressure, temperature, dewpoint, lapse_type='standard', params=None):
     r"""Calculate the profile a parcel takes through the atmosphere.
 
     The parcel starts at `temperature`, and `dewpoint`, lifted up
@@ -934,6 +1013,29 @@ def parcel_profile(pressure, temperature, dewpoint):
     dewpoint : `pint.Quantity`
         Starting dewpoint
 
+    lapse_type : `string`, optional
+        Definition of moist adiabat to use; if not given, it defaults to moist_lapse
+        Options:
+            'standard' for simplified pseudoadiabatic process
+            'pseudoadiabatic' for pseudoadiabatic moist process
+            'reversible' for reversible moist process
+            'so13' for Singh and O'Gorman (2013);  https://doi.org/10.1002/grl.50796
+            'r14' for Romps (2014); https://doi.org/10.1175/JCLI-D-14-00255.1
+        More info: https://glossary.ametsoc.org/wiki/Adiabatic_lapse_rate
+
+    params : `dict` or None, optional
+        External parameters used for the some lapse_types
+        Required parameters:
+            For 'so13': {
+                'ep0': entrainment constant [unitless],
+                'rh0': ambient relative humidity [unitless],
+                }
+            For 'r14': {
+                'de': scalar or 1-d array, detrainment rate [m**-1],
+                'ep': scalar or 1-d array, entrainment rate [m**-1],
+                'pa': 1-d array, optional, pressure levels defining detrainment and entrainment profile [Pa]
+                }
+
     Returns
     -------
     `pint.Quantity`
@@ -986,7 +1088,7 @@ def parcel_profile(pressure, temperature, dewpoint):
        Renamed ``dewpt`` parameter to ``dewpoint``
 
     """
-    _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint)
+    _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint, lapse_type, params)
     return concatenate((t_l, t_u))
 
 
@@ -1168,7 +1270,7 @@ def _check_pressure_error(pressure):
                                    'your sounding. Using scipy.signal.medfilt may fix this.')
 
 
-def _parcel_profile_helper(pressure, temperature, dewpoint):
+def _parcel_profile_helper(pressure, temperature, dewpoint, lapse_type, params):
     """Help calculate parcel profiles.
 
     Returns the temperature and pressure, above, below, and including the LCL. The
@@ -1205,7 +1307,7 @@ def _parcel_profile_helper(pressure, temperature, dewpoint):
                        'Output profile includes duplicate temperatures as a result.')
 
     # Find moist pseudo-adiabatic profile starting at the LCL, reversing above sorting
-    temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units)
+    temp_upper = moist_lapse(unique[::-1], temp_lower[-1], lapse_type=lapse_type, params=params).to(temp_lower.units)
     temp_upper = temp_upper[::-1][indices]
 
     # Return profile pieces

src/metpy/constants/nounit.py

@@ -7,8 +7,11 @@
 from ..units import units
 
 Rd = default.Rd.m_as('m**2 / K / s**2')
+Rv = default.Rv.m_as('m**2 / K / s**2')
 Lv = default.Lv.m_as('m**2 / s**2')
 Cp_d = default.Cp_d.m_as('m**2 / K / s**2')
+Cv_d = default.Cv_d.m_as('m**2 / K / s**2')
+Cp_l = default.Cp_l.m_as('m**2 / K / s**2')
 zero_degc = units.Quantity(0., 'degC').m_as('K')
 sat_pressure_0c = default.sat_pressure_0c.m_as('Pa')
 epsilon = default.epsilon.m_as('')