Skip to content

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)

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

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 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
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) -> 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.

    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_upper2 > p_upper):
        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

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 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
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