Lapse Integration (Simple)

fitting_2_layer(temp_env_lev, p_lev, temp_env_lower, p_lower, rh_lower, temp_env_upper, p_upper, temp_env_upper2, p_upper2, method_layer1='const', method_layer2='const', mod_parcel_method='add', n_lev_above_upper2_integral=0, temp_surf_lcl_calc=300, sanity_check=False, p_range_calc=2000)

Applies const_lapse_fitting or mod_parcel_lapse_fitting to each layer.

Different from lapse_integral.py in that it computes parcel temp, \(T_p(p)\) by equating \(MSE(T_{p,s}, r_s, p_s)\) rather than \(MSE^*(T_{LCL})\) to \(MSE^*(T_p(p))\). Where \(T_{p,s}\) is the parcel temperature at the surface starting from environmental LCL temperature.

Parameters:

Name	Type	Description	Default
temp_env_lev	Union[DataArray, ndarray]	[n_lev] Environmental temperature at pressure p_lev.	required
p_lev	Union[DataArray, ndarray]	[n_lev] Model pressure levels in Pa.	required
temp_env_lower	float	Environmental temperature at lower pressure level (nearer surface) p_lower.	required
p_lower	float	Pressure level to start profile (near surface).	required
rh_lower	float	Environmental relative humidity at p_lower.	required
temp_env_upper	float	Environmental temperature at the upper pressure level of the first layer (layer closest to surface) p_upper.	required
p_upper	float	Pressure level to end profile of first layer (layer closest to surface).	required
temp_env_upper2	float	Environmental temperature at the upper pressure level of the second layer (layer furthest from surface) p_upper2.	required
p_upper2	float	Pressure level to end profile of second layer (layer furthest from surface).	required
temp_parcel_lev		[n_lev] Parcel temperature (following moist adiabat) at pressure p_lev. Only required if either method_layer1 or method_layer2 are 'mod_parcel'.	required
temp_parcel_lower		Parcel temperature at lower pressure level (nearer surface) p_lower. Only required if either method_layer1 = 'mod_parcel'.	required
temp_parcel_upper		Parcel temperature at the upper pressure level of the first layer (layer closest to surface) p_upper. Only required if either method_layer1 or method_layer2 are 'mod_parcel'.	required
temp_parcel_upper2		Parcel temperature at the upper pressure level of the second layer (layer furthest from surface) p_upper2. Only required if either method_layer2 = 'mod_parcel'.	required
method_layer1	Literal['const']	Which fitting method to use for layer 1.	'const'
method_layer2	Literal['const', 'mod_parcel']	Which fitting method to use for layer 2.	'const'
n_lev_above_upper2_integral	int	Will return integral and integral_error not up to p_upper2 but up to the model level pressure n_lev_above_upper2_integral further from the surface than p_upper2. If n_lev_above_upper2_integral=0, upper limit of integral will be p_upper2.	0
temp_surf_lcl_calc	float	Surface temperature to use when computing \(\sigma_{LCL}\).	300
sanity_check	bool	If True will print a sanity check to ensure the calculation is correct.	False
p_range_calc	float	Only compute integral if p_upper2 > p_upper - p_range_calc. Units of Pa.	2000

Returns:

Name	Type	Description
lapse	ndarray	Lapse rate info for each layer. Bulk lapse rate if method_layer='const' or lapse rate adjustment if method_layer='mod_parcel'. Units are K/km.
integral	ndarray	Result of integral \(\int_{p_1}^{p_2} \Gamma_{env}(p) d\ln p\) of each layer. Units are K/km.
integral_error	ndarray	Result of integral \(\int_{p_1}^{p_2} |\Gamma_{env}(p) - \Gamma_{approx}| d\ln p\) of each layer. Units are K/km.

Source code in isca_tools/thesis/lapse_integral_simple.py

def fitting_2_layer(temp_env_lev: Union[xr.DataArray, np.ndarray], p_lev: Union[xr.DataArray, np.ndarray],
                    temp_env_lower: float, p_lower: float, rh_lower: float,
                    temp_env_upper: float, p_upper: float,
                    temp_env_upper2: float, p_upper2: float,
                    method_layer1: Literal['const'] = 'const',
                    method_layer2: Literal['const', 'mod_parcel'] = 'const',
                    mod_parcel_method: Literal['add', 'multiply'] = 'add',
                    n_lev_above_upper2_integral: int = 0,
                    temp_surf_lcl_calc: float = 300,
                    sanity_check: bool = False,
                    p_range_calc: float = 2000) -> Tuple[
    np.ndarray, np.ndarray, np.ndarray]:
    """
    Applies `const_lapse_fitting` or `mod_parcel_lapse_fitting` to each layer.

    Different from *lapse_integral.py* in that it computes parcel temp, $T_p(p)$ by equating
    $MSE(T_{p,s}, r_s, p_s)$ rather than $MSE^*(T_{LCL})$ to $MSE^*(T_p(p))$.
    Where $T_{p,s}$ is the parcel temperature at the surface starting from environmental LCL temperature.

    Args:
        temp_env_lev: `[n_lev]` Environmental temperature at pressure `p_lev`.
        p_lev: `[n_lev]` Model pressure levels in Pa.
        temp_env_lower: Environmental temperature at lower pressure level (nearer surface) `p_lower`.
        p_lower: Pressure level to start profile (near surface).
        rh_lower: Environmental relative humidity at `p_lower`.
        temp_env_upper: Environmental temperature at the upper pressure level of the first layer (layer closest to surface) `p_upper`.
        p_upper: Pressure level to end profile of first layer (layer closest to surface).
        temp_env_upper2: Environmental temperature at the upper pressure level of the second layer (layer furthest from surface) `p_upper2`.
        p_upper2: Pressure level to end profile of second layer (layer furthest from surface).
        temp_parcel_lev: `[n_lev]` Parcel temperature (following moist adiabat) at pressure `p_lev`.</br>
            Only required if either `method_layer1` or `method_layer2` are `'mod_parcel'`.
        temp_parcel_lower: Parcel temperature at lower pressure level (nearer surface) `p_lower`.</br>
            Only required if either `method_layer1 = 'mod_parcel'`.
        temp_parcel_upper: Parcel temperature at the upper pressure level of the first layer (layer closest to surface) `p_upper`.</br>
            Only required if either `method_layer1` or `method_layer2` are `'mod_parcel'`.
        temp_parcel_upper2: Parcel temperature at the upper pressure level of the second layer (layer furthest from surface) `p_upper2`.</br>
            Only required if either `method_layer2 = 'mod_parcel'`.
        method_layer1: Which fitting method to use for layer 1.
        method_layer2: Which fitting method to use for layer 2.
        n_lev_above_upper2_integral: Will return `integral` and `integral_error` not up to `p_upper2` but up to
            the model level pressure `n_lev_above_upper2_integral` further from the surface than `p_upper2`.
            If `n_lev_above_upper2_integral=0`, upper limit of integral will be `p_upper2`.
        temp_surf_lcl_calc: Surface temperature to use when computing $\sigma_{LCL}$.
        sanity_check: If `True` will print a sanity check to ensure the calculation is correct.
        p_range_calc: Only compute integral if `p_upper2 > p_upper - p_range_calc`. Units of Pa.

    Returns:
        lapse: Lapse rate info for each layer. Bulk lapse rate if `method_layer='const'` or lapse rate adjustment
            if `method_layer='mod_parcel'`. Units are *K/km*.
        integral: Result of integral $\int_{p_1}^{p_2} \Gamma_{env}(p) d\ln p$ of each layer. Units are *K/km*.
        integral_error: Result of integral $\int_{p_1}^{p_2} |\Gamma_{env}(p) - \Gamma_{approx}| d\ln p$ of each layer.
            Units are *K/km*.
    """
    if (p_upper > p_lower-p_range_calc) | (p_upper2 > p_upper-p_range_calc):
        # LCL or Surface too close to FT level to compute
        return np.asarray([np.nan, np.nan]), np.asarray([np.nan, np.nan]), np.asarray([np.nan, np.nan])
    if method_layer1 == 'const':
        lapse1, lapse_integral1, lapse_integral_error1 = \
            const_lapse_fitting(temp_env_lev, p_lev, temp_env_lower, p_lower, temp_env_upper, p_upper,
                                sanity_check=sanity_check)
    else:
        raise ValueError(f'method_layer1 = {method_layer1} not recognized.')

    if method_layer2 == 'const':
        lapse2, lapse_integral2, lapse_integral_error2 = \
            const_lapse_fitting(temp_env_lev, p_lev, temp_env_upper, p_upper, temp_env_upper2, p_upper2,
                                n_lev_above_upper2_integral, sanity_check=sanity_check)
    elif method_layer2 == 'mod_parcel':
        # Compute parcel temp at p_lower following a dry adiabat from p_upper
        temp_parcel_lower = dry_profile_temp(temp_env_upper, p_upper, p_lower)
        lapse2, lapse_integral2, lapse_integral_error2 = \
            mod_parcel_lapse_fitting(temp_env_lev, p_lev, temp_env_upper, p_upper, temp_env_upper2, p_upper2,
                                     temp_parcel_surf=temp_parcel_lower, rh_parcel_surf=rh_lower, p_surf=p_lower,
                                     n_lev_above_upper_integral=n_lev_above_upper2_integral,
                                     temp_surf_lcl_calc=temp_surf_lcl_calc,
                                     method=mod_parcel_method, sanity_check=sanity_check)
    else:
        raise ValueError(f'method_layer2 = {method_layer2} not recognized.')

    lapse = np.asarray([lapse1, lapse2])
    lapse_integral = np.asarray([lapse_integral1, lapse_integral2])
    lapse_integral_error = np.asarray([lapse_integral_error1, lapse_integral_error2])
    return lapse, lapse_integral, lapse_integral_error

get_lnb_ind(temp_env_lev, p_lev, rh_surf, lapse_mod_D=0, p_surf=None, temp_surf=None, temp_surf_lcl_calc=300)

Returns a proxy for the LNB (level of neutral buoyancy): the index closest to space where T_parcel - T_env is positive i.e. buoyant.

Parameters:

Name	Type	Description	Default
temp_env_lev	ndarray	float [n_lev] Environmental temperature on model levels. Units: Kelvin.	required
p_lev	ndarray	float [n_lev] Pressure on model levels. Units: Pa.	required
rh_surf	float	Relative humidity at p_surf, used to calculate LCL.	required
lapse_mod_D	float	The parameter, such that boundary layer lapse rate equals lapse_dry + lapse_mod_D. Units: K/m. If zero more classical definition of \(T_{parc}\) used i.e. parcel rising from environmental surface temperature. If provide modification to give environmental lapse rate, a \(T_{parc-L}\) version used with parcel rising from environmental LCL temperature.	0
p_surf	Optional[float]	Surface pressure used to calculate LCL. Units: Pa. If not provided, will set to model level closest to the surface.	None
temp_surf	Optional[float]	Temperature at p_surf. Units: K. If not provided, will set to model level closest to the surface.	None
temp_surf_lcl_calc		Surface temperature to use when computing \(\sigma_{LCL}\).	300

Returns:

Type	Description
int	Index of model level closest to the top of the atmosphere for which parcel temperature is warmer than environmental
int	temperature. I.e. index of the buoyant layer furthest from the surface.
int	If no buoyant layer, it will return the index of surface.

Source code in isca_tools/thesis/lapse_integral_simple.py

def get_lnb_ind(temp_env_lev: np.ndarray, p_lev: np.ndarray, rh_surf: float, lapse_mod_D: float=0,
                p_surf: Optional[float]=None,
                temp_surf:Optional[float]=None,
                temp_surf_lcl_calc=300)-> int:
    """
    Returns a proxy for the LNB (level of neutral buoyancy): the index closest to space where
    `T_parcel - T_env` is positive i.e. buoyant.

    Args:
        temp_env_lev: `float [n_lev]`</br>
            Environmental temperature on model levels. Units: Kelvin.
        p_lev: `float [n_lev]`</br>
            Pressure on model levels. Units: Pa.
        rh_surf: Relative humidity at `p_surf`, used to calculate LCL.
        lapse_mod_D: The parameter, such that boundary layer lapse rate equals `lapse_dry + lapse_mod_D`.
            Units: K/m.</br>
            If zero more classical definition of $T_{parc}$ used i.e. parcel rising from environmental surface temperature.
            If provide modification to give environmental lapse rate, a $T_{parc-L}$ version used with
            parcel rising from environmental LCL temperature.
        p_surf: Surface pressure used to calculate LCL. Units: Pa.
            If not provided, will set to model level closest to the surface.
        temp_surf: Temperature at `p_surf`. Units: K.
            If not provided, will set to model level closest to the surface.
        temp_surf_lcl_calc: Surface temperature to use when computing $\sigma_{LCL}$.

    Returns:
        Index of model level closest to the top of the atmosphere for which parcel temperature is warmer than environmental
        temperature. I.e. index of the buoyant layer furthest from the surface.
        If no buoyant layer, it will return the index of surface.
    """
    surf_ind = np.argmax(p_lev)
    if surf_ind == 0:
        raise ValueError('Must have pressure ascending in p_lev')
    if p_surf is None:
        p_surf = p_lev[surf_ind]
        if temp_surf is not None:
            raise ValueError('If provide p_surf, must also provide temp_surf.')
        temp_surf = temp_env_lev[surf_ind]
    elif temp_surf is not None:
        raise ValueError('If do not provide p_surf, must also not provide temp_surf.')
    p_lcl = lcl_sigma_bolton_simple(rh_surf, temp_surf_lcl_calc) * p_surf
    # Only Consider pressure values above LCL
    mask_ind = np.where(p_lev<p_lcl)[0]
    temp_env_lev = temp_env_lev[p_lev<p_lcl]
    p_lev = p_lev[p_lev<p_lcl]
    lnb_ind = surf_ind    # default value of lnb_ind if always buoyant i.e. in space
    for i in range(temp_env_lev.size):
        temp_parcel = get_temp_mod_parcel(rh_surf, p_surf, p_lev[i], lapse_mod_D,
                                          0, temp_surf, lapse_coords='z')
        if temp_parcel > temp_env_lev[i]:
            lnb_ind = mask_ind[i]
            break
    return lnb_ind

mod_parcel_lapse_fitting(temp_env_lev, p_lev, temp_env_lower, p_lower, temp_env_upper, p_upper, temp_parcel_surf, rh_parcel_surf, p_surf, n_lev_above_upper_integral=0, temp_surf_lcl_calc=300, sanity_check=False, method='add')

Find the constant, \(\eta\) that needs adding to parcel lapse rate such that \(\int_{p_1}^{p_2} \Gamma_{env}(p) d\ln p = \int_{p_1}^{p_2} \Gamma_p(p, T_p(p)) + \eta d\ln p\). Then computes the error in this approximation: \(\int_{p_1}^{p_2} |\Gamma_{env}(p) - \Gamma_p(p, T_p(p)) - \eta| d\ln p\).

where \(\Gamma_p(p, T)\) is the parcel (moist adiabatic) lapse rate and \(T_p(p)\) is the parcel temperature at pressure \(p\) starting at \(p_1\).

Different from lapse_integral.py in that it computes parcel temp, \(T_p(p)\) by equating \(MSE(T_{p,s}, r_s, p_s)\) rather than \(MSE^*(T_{LCL})\) to \(MSE^*(T_p(p))\). Where \(T_{p,s}\) is the parcel temperature at the surface starting from environmental LCL temperature.

Parameters:

Name	Type	Description	Default
temp_env_lev	ndarray	[n_lev] Environmental temperature at pressure p_lev.	required
p_lev	ndarray	[n_lev] Model pressure levels in Pa.	required
temp_env_lower	float	Environmental temperature at lower pressure level (nearer surface) p_lower.	required
p_lower	float	Pressure level to start profile (near surface).	required
temp_env_upper	float	Environmental temperature at upper pressure level (further from surface) p_upper.	required
p_upper	float	Pressure level to end profile (further from surface).	required
temp_parcel_surf	float	Temperature of the parcel at p_surf.	required
rh_parcel_surf	float	Relative humidity of the parcel at p_surf.	required
p_surf	float	Pressure at p_surf. Units: Pa.	required
n_lev_above_upper_integral	int	Will return integral and integral_error not up to p_upper but up to the model level pressure n_lev_above_upper_integral further from the surface than p_upper. If n_lev_above_upper_integral=0, upper limit of integral will be p_upper.	0
temp_surf_lcl_calc	float	Surface temperature to use when computing \(\sigma_{LCL}\).	300
sanity_check	bool	If True will print a sanity check to ensure the calculation is correct.	False

Returns:

Name	Type	Description
lapse_diff_const	float	Lapse rate adjustment, \(\eta\) which needs to be added to \(\Gamma_{p}(p, T_p(p))\) so integral matches that of environmental lapse rate. Units are K/km.
integral	float	Result of integral \(\int_{p_1}^{p_2} \Gamma_{env}(p) d\ln p\). Units are K/km.
integral_error	float	Result of integral \(\int_{p_1}^{p_2} |\Gamma_{env}(p) - \Gamma_p(p) - \eta| d\ln p\). Units are K/km.

Source code in isca_tools/thesis/lapse_integral_simple.py

def mod_parcel_lapse_fitting(temp_env_lev: np.ndarray,
                             p_lev: np.ndarray,
                             temp_env_lower: float, p_lower: float,
                             temp_env_upper: float, p_upper: float,
                             temp_parcel_surf: float,
                             rh_parcel_surf: float,
                             p_surf: float,
                             n_lev_above_upper_integral: int = 0,
                             temp_surf_lcl_calc: float = 300,
                             sanity_check: bool = False,
                             method: Literal['add', 'multiply'] = 'add') -> Tuple[
    float, float, float]:
    """
    Find the constant, $\eta$ that needs adding to parcel lapse rate such that
    $\int_{p_1}^{p_2} \Gamma_{env}(p) d\ln p = \int_{p_1}^{p_2} \Gamma_p(p, T_p(p)) + \eta d\ln p$.
    Then computes the error in this approximation:
    $\int_{p_1}^{p_2} |\Gamma_{env}(p) - \Gamma_p(p, T_p(p)) - \eta| d\ln p$.

    where $\Gamma_p(p, T)$ is the parcel (moist adiabatic) lapse rate and $T_p(p)$ is the parcel temperature
    at pressure $p$ starting at $p_1$.

    Different from *lapse_integral.py* in that it computes parcel temp, $T_p(p)$ by equating
    $MSE(T_{p,s}, r_s, p_s)$ rather than $MSE^*(T_{LCL})$ to $MSE^*(T_p(p))$.
    Where $T_{p,s}$ is the parcel temperature at the surface starting from environmental LCL temperature.

    Args:
        temp_env_lev: `[n_lev]` Environmental temperature at pressure `p_lev`.
        p_lev: `[n_lev]` Model pressure levels in Pa.
        temp_env_lower: Environmental temperature at lower pressure level (nearer surface) `p_lower`.
        p_lower: Pressure level to start profile (near surface).
        temp_env_upper: Environmental temperature at upper pressure level (further from surface) `p_upper`.
        p_upper: Pressure level to end profile (further from surface).
        temp_parcel_surf: Temperature of the parcel at `p_surf`.
        rh_parcel_surf: Relative humidity of the parcel at `p_surf`.
        p_surf: Pressure at `p_surf`. Units: Pa.
        n_lev_above_upper_integral: Will return `integral` and `integral_error` not up to `p_upper` but up to
            the model level pressure `n_lev_above_upper_integral` further from the surface than `p_upper`.
            If `n_lev_above_upper_integral=0`, upper limit of integral will be `p_upper`.
        temp_surf_lcl_calc: Surface temperature to use when computing $\sigma_{LCL}$.
        sanity_check: If `True` will print a sanity check to ensure the calculation is correct.

    Returns:
        lapse_diff_const: Lapse rate adjustment, $\eta$ which needs to be added to $\Gamma_{p}(p, T_p(p))$
            so integral matches that of environmental lapse rate. Units are *K/km*.
        integral: Result of integral $\int_{p_1}^{p_2} \Gamma_{env}(p) d\ln p$. Units are *K/km*.
        integral_error: Result of integral $\int_{p_1}^{p_2} |\Gamma_{env}(p) - \Gamma_p(p) - \eta| d\ln p$.
            Units are *K/km*.
    """
    if n_lev_above_upper_integral == 0:
        p_integ_upper = p_upper
    else:
        if n_lev_above_upper_integral < 0:
            raise ValueError('n_lev_above_upper_integral must be greater than 0')
        if p_lev[1] < p_lev[0]:
            raise ValueError('p_lev[1] must be greater than p_lev[0]')
        ind_upper = np.where(p_lev < p_upper)[0][-n_lev_above_upper_integral]
        p_integ_upper = p_lev[ind_upper]

    # Get parcel upper temp following moist adiabat, using surface MSE
    temp_parcel_upper = get_temp_mod_parcel(rh_parcel_surf, p_surf, p_upper, 0, 0,
                                            temp_parcel_surf, temp_surf_lcl_calc=temp_surf_lcl_calc)

    if method == 'multiply':
        lapse_diff_const = np.log(temp_env_upper / temp_env_lower) / np.log(temp_parcel_upper / temp_env_lower) - 1
    elif method == 'add':
        # Compute integral of deviation between environmental and parcel lapse rate between p_lower and p_upper
        lapse_diff_integral = g / R * (
                np.log(temp_env_upper / temp_env_lower) - np.log(temp_parcel_upper / temp_env_lower))
        # Compute the constant needed to be added to the parcel lapse rate at each level to make above integral equal 0.
        lapse_diff_const = lapse_diff_integral / np.log(p_upper / p_lower)

    temp_env_approx_lev = np.zeros_like(temp_env_lev)
    for i in range(temp_env_approx_lev.size):
        if (p_lev[i] > p_lower) | (p_lev[i] < p_integ_upper):
            # If not in the pressure range required for integral, then keep as zero
            continue
        temp_env_approx_lev[i] = get_temp_mod_parcel(rh_parcel_surf, p_surf, p_lev[i], 0, lapse_diff_const,
                                                     temp_parcel_surf, temp_surf_lcl_calc=temp_surf_lcl_calc,
                                                     method=method)

    if sanity_check:
        temp_env_approx_upper = get_temp_mod_parcel(rh_parcel_surf, p_surf, p_upper, 0, lapse_diff_const,
                                                    temp_parcel_surf, temp_surf_lcl_calc=temp_surf_lcl_calc,
                                                    method=method)
        lapse_integral = g / R * np.log(temp_env_upper / temp_env_lower)
        # sanity check, this should be the same as lapse_integral
        lapse_integral_approx = integral_lapse_dlnp_hydrostatic(temp_env_approx_lev, p_lev, p_lower, p_upper,
                                                                temp_env_lower, temp_env_upper)
        print(
            f'Actual lapse integral: {lapse_integral * 1000:.3f} K/km\nApprox lapse integral: {lapse_integral_approx * 1000:.3f} K/km')
        # Will use lapse rate such that approx value of T_upper is exact. Check that here
        print(f'Actual temp_upper: {temp_env_upper:.3f} K\nApprox temp_upper: {temp_env_approx_upper:.3f} K')

    # Compute error in the integral
    if n_lev_above_upper_integral == 0:
        lapse_integral, lapse_integral_error = integral_and_error_calc(temp_env_lev, temp_env_approx_lev, p_lev,
                                                                       temp_env_lower, temp_env_lower, p_lower,
                                                                       temp_env_upper, temp_env_upper, p_upper)
    else:
        lapse_integral, lapse_integral_error = integral_and_error_calc(temp_env_lev, temp_env_approx_lev, p_lev,
                                                                       temp_env_lower, temp_env_lower, p_lower,
                                                                       float(temp_env_lev[ind_upper]),
                                                                       float(temp_env_approx_lev[ind_upper]),
                                                                       float(p_lev[ind_upper]))
    return lapse_diff_const * (1000 if method == 'add' else 1), lapse_integral * 1000, lapse_integral_error * 1000