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', const_layer1_method='bulk', mod_parcel_method='add', force_parcel=False, 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
method_layer1	Literal['const']	Which fitting method to use for layer 1. Not required if force_parcel=True.	'const'
method_layer2	Literal['const', 'mod_parcel']	Which fitting method to use for layer 2. Not required if force_parcel=True.	'const'
const_layer1_method	Literal['bulk', 'optimal']	If 'bulk', will compute lapse rate which intersects exactly with temp_env_lower at p_lower and temp_env_upper at p_upper. If optimal, will use fit_lapse_L1_dlnT to compute the optimal lapse rate i.e., to minimize the error. Will still intersect temp_env_lower at p_lower, but not temp_env_upper at p_upper.	'bulk'
mod_parcel_method	Literal['add', 'multiply']	Recommend add, multiply is another method I tried where you don't add to the moist adiabatic lapse rate at each level but multiply it by a scaling factor.	'add'
force_parcel	bool	If True, will use parcel_fitting to find error and integral values. I.e. force dry adiabat below surface and moist adiabat above. But parcel rising from surface so LCL parcel temperature will differ from environmental.	False
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',
                    const_layer1_method: Literal['bulk', 'optimal'] = 'bulk',
                    mod_parcel_method: Literal['add', 'multiply'] = 'add',
                    force_parcel: bool = False,
                    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).
        method_layer1: Which fitting method to use for layer 1. Not required if `force_parcel=True`.
        method_layer2: Which fitting method to use for layer 2. Not required if `force_parcel=True`.
        const_layer1_method: If 'bulk', will compute lapse rate which intersects exactly with
            `temp_env_lower` at `p_lower` and `temp_env_upper` at `p_upper`.
            If `optimal`, will use `fit_lapse_L1_dlnT` to compute the optimal lapse rate i.e., to minimize the error.
            Will still intersect `temp_env_lower` at `p_lower`, but not `temp_env_upper` at `p_upper`.
        mod_parcel_method: Recommend `add`, `multiply` is another method I tried where you don't add to the moist
            adiabatic lapse rate at each level but multiply it by a scaling factor.
        force_parcel: If `True`, will use `parcel_fitting` to find error and integral values.
            I.e. force dry adiabat below surface and moist adiabat above. But parcel rising from surface
            so LCL parcel temperature will differ from environmental.
        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) | (rh_lower <= 0) | (rh_lower >= 1):
        # (LCL too close to surface) or (LCL 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 force_parcel:
        return parcel_fitting(temp_env_lev, p_lev, temp_env_lower, p_lower, rh_lower, temp_env_upper, p_upper,
                              temp_env_upper2, p_upper2, n_lev_above_upper2_integral, temp_surf_lcl_calc)
    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, lapse_method=const_layer1_method)
    else:
        raise ValueError(f'method_layer1 = {method_layer1} not recognized.')

    if const_layer1_method == 'optimal':
        temp_approx_upper = get_temp_const_lapse(p_upper, temp_env_lower, p_lower, lapse1 / 1000)
    else:
        temp_approx_upper = temp_env_upper

    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,
                                lapse_method='bulk', temp_approx_lower=temp_approx_upper)
    elif method_layer2 == 'mod_parcel':
        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_lower=temp_approx_upper, 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.
int	Will be -1 if could not be computed.

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.
        Will be -1 if could not be computed.
    """
    if (rh_surf <= 0) or (rh_surf >= 1) or np.isnan(lapse_mod_D):
        return -1
    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_lower, 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_lower	Optional[float]	Temperature of the parcel at p_lower. Will set to temp_env_lower if not provided. This is the temperature of the approximate profile at p_lower. Does not have to match temp_env_lower.	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
method	Literal['add', 'multiply']	Recommend add, multiply is another method I tried where you don't add to the moist adiabatic lapse rate at each level but multiply it by a scaling factor.	'add'

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_lower: Optional[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_lower: Temperature of the parcel at `p_lower`. Will set to `temp_env_lower` if not provided.
            This is the temperature of the approximate profile at `p_lower`. Does not have to match `temp_env_lower`.
        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.
        method: Recommend `add`, `multiply` is another method I tried where you don't add to the moist
            adiabatic lapse rate at each level but multiply it by a scaling factor.

    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 temp_parcel_lower is None:
        temp_parcel_lower = temp_env_lower
    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]
    # Compute parcel temp at surface following a dry adiabat from the parcel temp at p_lower
    temp_parcel_surf = dry_profile_temp(temp_parcel_lower, p_lower, p_surf)

    # 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_parcel_lower) / np.log(
            temp_parcel_upper / temp_parcel_lower) - 1
    elif method == 'add':
        # Note temp_parcel_lower cancels out so not needed to compute the value of lapse_diff_const. Just include it for clarity
        # 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_parcel_lower) - np.log(temp_parcel_upper / temp_parcel_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_approx_lev = np.zeros_like(temp_env_lev)
    for i in range(temp_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_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_parcel_lower)
        # sanity check, this should be the same as lapse_integral
        lapse_integral_approx = integral_lapse_dlnp_hydrostatic(temp_approx_lev, p_lev, p_lower, p_upper,
                                                                temp_parcel_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:
        # By definition temp_approx = temp_env at p_upper but at p_lower temp_approx=temp_parcel_lower
        lapse_integral, lapse_integral_error = integral_and_error_calc(temp_env_lev, temp_approx_lev, p_lev,
                                                                       temp_env_lower, temp_parcel_lower, p_lower,
                                                                       temp_env_upper, temp_env_upper, p_upper)
    else:
        # Note that at p_lower temp_approx=temp_parcel_lower
        lapse_integral, lapse_integral_error = integral_and_error_calc(temp_env_lev, temp_approx_lev, p_lev,
                                                                       temp_env_lower, temp_parcel_lower, p_lower,
                                                                       float(temp_env_lev[ind_upper]),
                                                                       float(temp_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

parcel_fitting(temp_env_lev, p_lev, temp_env_lower, p_lower, rh_lower, temp_env_upper, p_upper, temp_env_upper2, p_upper2, n_lev_above_upper2_integral=0, temp_surf_lcl_calc=300, return_temp_parcel_lev=False)

If follow parcel rising from the surface, this returns the error and integral value, of this profile relative to the environmental profile.

Parameters:

Name	Type	Description	Default
temp_env_lev	ndarray		required
p_lev	ndarray	Surface pressure	required
temp_env_lower	float		required
p_lower	float		required
rh_lower	float		required
temp_env_upper	float		required
p_upper	float	LCL pressure	required
temp_env_upper2	float		required
p_upper2	float	FT pressure	required
n_lev_above_upper2_integral	int		0
temp_surf_lcl_calc	float		300

Returns:

Name	Type	Description
lapse		[lapse_dry, 0] i.e. dry adiabat below LCL and moist adiabat above. 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.

Source code in isca_tools/thesis/lapse_integral_simple.py

def parcel_fitting(temp_env_lev: np.ndarray, p_lev: 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,
                   n_lev_above_upper2_integral: int = 0,
                   temp_surf_lcl_calc: float = 300,
                   return_temp_parcel_lev: bool = False):
    """
    If follow parcel rising from the surface, this returns the error and integral value, of this profile
    relative to the environmental profile.

    Args:
        temp_env_lev:
        p_lev: Surface pressure
        temp_env_lower:
        p_lower:
        rh_lower:
        temp_env_upper:
        p_upper: LCL pressure
        temp_env_upper2:
        p_upper2: FT pressure
        n_lev_above_upper2_integral:
        temp_surf_lcl_calc:

    Returns:
        lapse: `[lapse_dry, 0]` i.e. dry adiabat below LCL and moist adiabat above. 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*.
    """
    temp_parcel_upper = dry_profile_temp(temp_env_lower, p_lower, p_upper)  # parcel temp at lcl
    temp_parcel_lev = np.zeros_like(temp_env_lev)
    for i in range(p_lev.size):
        if p_lev[i] >= p_upper:
            # Dry adiabatic region - just follow dry adiabat
            temp_parcel_lev[i] = dry_profile_temp(temp_env_lower, p_lower, p_lev[i])
        else:
            # Moist adiabatic region - same modMSE as surface
            temp_parcel_lev[i] = get_temp_mod_parcel(rh_lower, p_lower, p_lev[i], 0, 0,
                                                     temp_env_lower, temp_surf_lcl_calc=temp_surf_lcl_calc)
    lapse_integral1, lapse_integral_error1 = integral_and_error_calc(temp_env_lev, temp_parcel_lev, p_lev,
                                                                     temp_env_lower, temp_env_lower, p_lower,
                                                                     temp_env_upper, temp_parcel_upper, p_upper)
    if n_lev_above_upper2_integral == 0:
        temp_parcel_upper2 = get_temp_mod_parcel(rh_lower, p_lower, p_upper2, 0, 0,
                                                 temp_env_lower, temp_surf_lcl_calc=temp_surf_lcl_calc)
        lapse_integral2, lapse_integral_error2 = integral_and_error_calc(temp_env_lev, temp_parcel_lev, p_lev,
                                                                         temp_env_upper, temp_parcel_upper, p_upper,
                                                                         temp_env_upper2, temp_parcel_upper2, p_upper2)
    else:
        ind_upper = np.where(p_lev < p_upper2)[0][-n_lev_above_upper2_integral]
        lapse_integral2, lapse_integral_error2 = integral_and_error_calc(temp_env_lev, temp_parcel_lev, p_lev,
                                                                         temp_env_upper, temp_parcel_upper, p_upper,
                                                                         temp_env_lev[ind_upper],
                                                                         temp_parcel_lev[ind_upper],
                                                                         p_lev[ind_upper])
    lapse_diff_const = np.asarray([lapse_dry, 0])
    lapse_integral = np.asarray([lapse_integral1, lapse_integral2])
    lapse_integral_error = np.asarray([lapse_integral_error1, lapse_integral_error2])
    if return_temp_parcel_lev:
        return lapse_diff_const * 1000, lapse_integral * 1000, lapse_integral_error * 1000, temp_parcel_lev
    else:
        return lapse_diff_const * 1000, lapse_integral * 1000, lapse_integral_error * 1000