Simple Betts-Miller

`get_temp_ref(temp_start, p_start, sphum_start, p_full, p_ref=100000.0)`

This replicates the way the reference temperature profile is computed in Isca with the Simple Betts-Miller convection scheme.

Below the LCL, it is set to a dry adiabat, and above it the ref_temp_above_lcl function is used, which is very similar to the moist adiabat.

As well as the reference temperature profile, the LCL temperature and pressure are also returned.

Parameters:

Name	Type	Description	Default
`temp_start`	`float`	Starting temperature of parcel, \(T_{start}\). Units: Kelvin.	required
`p_start`	`float`	Pressure, \(p_{start}\), in Pa corresponding to starting point of dry ascent i.e. near surface pressure.	required
`sphum_start`	`float`	Starting specific humidity of parcel, \(q_{start}\). Units: kg/kg.	required
`p_full`	`ndarray`	`float [n_p_levels]`. Full model pressure levels in ascending order. `p_full[0]` represents space and `p_full[-1]` is the surface. Units: Pa.	required
`p_ref`	`float`	Reference pressure, \(p_{ref}\) in Pa. It is a parameter in the `qe_moist_convection` namelist.	`100000.0`

Returns:

Type	Description
`ndarray`	`temp_ref`: `float [n_p_levels]` The reference temperature at each pressure level in `p_full`.
`float`	`temp_lcl`: Temperature of LCL in K.
`float`	`p_lcl`: Pressure of LCL in Pa.

Source code in isca_tools/convection/simple_betts_miller.py

def get_temp_ref(temp_start: float, p_start: float, sphum_start: float,
                 p_full: np.ndarray, p_ref: float = 1e5) -> Tuple[np.ndarray, float, float]:
    """
    This replicates the way the reference temperature profile is computed in Isca with the
    [Simple Betts-Miller](https://jduffield65.github.io/Isca/namelists/convection/qe_moist_convection/)
    convection scheme.

    Below the LCL, it is set to a dry adiabat, and above it the `ref_temp_above_lcl` function is used, which is very
    similar to the moist adiabat.

    As well as the reference temperature profile, the LCL temperature and pressure are also returned.

    Args:
        temp_start: Starting temperature of parcel, $T_{start}$. Units: *Kelvin*.
        p_start: Pressure, $p_{start}$, in *Pa* corresponding to starting point of dry ascent i.e. near surface
            pressure.
        sphum_start: Starting specific humidity of parcel, $q_{start}$. Units: *kg/kg*.
        p_full: `float [n_p_levels]`.</br>
            Full model pressure levels in ascending order. `p_full[0]` represents space and `p_full[-1]` is
            the surface.</br>Units: *Pa*.
        p_ref: Reference pressure, $p_{ref}$ in *Pa*. It is a parameter in the `qe_moist_convection` namelist.

    Returns:
        `temp_ref`: `float [n_p_levels]`</br>
            The reference temperature at each pressure level in `p_full`.
        `temp_lcl`: Temperature of *LCL* in *K*.
        `p_lcl`: Pressure of *LCL* in *Pa*.
    """
    # Important that all variables are of correct datatype
    if not isinstance(temp_start, float):
        warnings.warn('Changing temp_start to a float')
        temp_start = float(temp_start)
    if not isinstance(p_start, float):
        warnings.warn('Changing p_start to a float')
        p_start = float(p_start)
    if not isinstance(sphum_start, float):
        warnings.warn('Changing sphum_start to a float')
        sphum_start = float(sphum_start)
    if not isinstance(p_ref, float):
        warnings.warn('Changing p_ref to a float')
        p_ref = float(p_ref)
    if not isinstance(p_full, np.ndarray):
        warnings.warn('Changing p_full to a numpy array')
        p_full = np.asarray(p_full)

    temp_lcl = lcl_temp(temp_start, p_start, sphum_start, p_ref)
    p_lcl = dry_profile_pressure(temp_start, p_start, temp_lcl)
    temp_ref = np.zeros(len(p_full))
    temp_ref[p_full >= p_lcl] = dry_profile_temp(temp_start, p_start, p_full[p_full >= p_lcl])
    temp_ref[p_full < p_lcl] = ref_temp_above_lcl(temp_lcl, p_lcl, p_full[p_full < p_lcl])
    return temp_ref, temp_lcl, p_lcl

`lcl_temp(temp_start, p_start, sphum_start, p_ref=100000.0)`

Function to replicate the way LCL temperature is computed in Isca with the Simple Betts-Miller convection scheme.

It is based on two properties of dry ascent:

Potential temperature is conserved so surface potential temperature = potential temperature at the \(LCL\):

\[\theta = \theta_{start}(T_{start}, p_{start}) = T_{start}\bigg(\frac{p_{ref}}{p_{start}}\bigg)^{\kappa} = \theta_{LCL}(T_{LCL}, p_{LCL}) = T_{LCL}\bigg(\frac{p_{ref}}{p_{LCL}}\bigg)^{\kappa}\]

Mixing ratio, \(w\), is conserved in unsaturated adiabatic ascent because there is no precipitation, and at the \(LCL\), \(w_{LCL} = w_{sat}\) because by definition at the \(LCL\), the air is saturated:

\[w = w_{start} = \frac{q_{start}}{1-q_{start}} = w_{sat} = \frac{\epsilon e_s(T_{LCL})}{p_{LCL}-e_s(T_{LCL})}\]

\(q\) is specific humidity, \(\epsilon = R_{dry}/R_v\) is the ratio of gas constant for dry air to vapour and \(\kappa = R_{dry}/c_p\). \(p_{ref}\) is taken to be \(100,000 Pa\) to be consistent with the value used in Isca.

So we have two equations for two unknowns, \(T_{LCL}\) and \(p_{LCL}\). By eliminating \(p_{LCL}\), we can get an equation where RHS is just a function of \(T_{LCL}\) and the LHS consists only of known quantities:

\[\theta(T_{start})^{\kappa}p_{ref} \frac{w(q_{start})}{w(q_{start}) + \epsilon} = T_{LCL}^{-1/\kappa}e_s(T_{LCL})\]

This can then be solved using Newton iteration to get the value of \(T_{LCL}\).

Parameters:

Name	Type	Description	Default
`temp_start`	`float`	Starting temperature of parcel, \(T_{start}\). Units: Kelvin.	required
`p_start`	`float`	Pressure, \(p_{start}\), in Pa corresponding to starting point of dry ascent i.e. near surface pressure.	required
`sphum_start`	`float`	Starting specific humidity of parcel, \(q_{start}\). Units: kg/kg.	required
`p_ref`	`float`	Reference pressure, \(p_{ref}\) in Pa. It is a parameter in the `qe_moist_convection` namelist.	`100000.0`

Returns:

Type	Description
`float`	Temperature of LCL in K.

Source code in isca_tools/convection/simple_betts_miller.py

def lcl_temp(temp_start: float, p_start: float, sphum_start: float, p_ref: float = 1e5) -> float:
    """
    Function to replicate the way
    [LCL temperature is computed](https://github.com/ExeClim/Isca/blob/9560521e1ba5ce27a13786ffdcb16578d0bd00da/src/
    atmos_param/qe_moist_convection/qe_moist_convection.F90#L1092-L1130)
    in *Isca* with the
    [Simple Betts-Miller](https://jduffield65.github.io/Isca/namelists/convection/qe_moist_convection/)
    convection scheme.

    It is based on two properties of dry ascent:

    Potential temperature is conserved so surface potential temperature = potential temperature at the $LCL$:

    $$\\theta = \\theta_{start}(T_{start}, p_{start}) = T_{start}\\bigg(\\frac{p_{ref}}{p_{start}}\\bigg)^{\kappa} =
    \\theta_{LCL}(T_{LCL}, p_{LCL}) = T_{LCL}\\bigg(\\frac{p_{ref}}{p_{LCL}}\\bigg)^{\kappa}$$

    Mixing ratio, $w$, is conserved in unsaturated adiabatic ascent because there is no precipitation,
    and at the $LCL$, $w_{LCL} = w_{sat}$ because by definition at the $LCL$, the air is saturated:

    $$w = w_{start} = \\frac{q_{start}}{1-q_{start}} = w_{sat} = \\frac{\epsilon e_s(T_{LCL})}{p_{LCL}-e_s(T_{LCL})}$$

    $q$ is specific humidity, $\epsilon = R_{dry}/R_v$ is the ratio of gas constant for dry air to vapour and
    $\kappa = R_{dry}/c_p$. $p_{ref}$ is taken to be $100,000 Pa$ to be consistent with the
    [value used in Isca](https://github.com/ExeClim/Isca/blob/9560521e1ba5ce27a13786ffdcb16578d0bd00da/src/
    atmos_param/qe_moist_convection/qe_moist_convection.F90#L74-L75).

    So we have two equations for two unknowns, $T_{LCL}$ and $p_{LCL}$. By eliminating $p_{LCL}$,
    we can get an equation where RHS is just a function of $T_{LCL}$ and the LHS consists only of known quantities:

    $$\\theta(T_{start})^{\kappa}p_{ref} \\frac{w(q_{start})}{w(q_{start}) + \epsilon} =
    T_{LCL}^{-1/\kappa}e_s(T_{LCL})$$

    This can then be solved using Newton iteration to get the value of $T_{LCL}$.

    Args:
        temp_start: Starting temperature of parcel, $T_{start}$. Units: *Kelvin*.
        p_start: Pressure, $p_{start}$, in *Pa* corresponding to starting point of dry ascent i.e. near surface
            pressure.
        sphum_start: Starting specific humidity of parcel, $q_{start}$. Units: *kg/kg*.
        p_ref: Reference pressure, $p_{ref}$ in *Pa*. It is a parameter in the `qe_moist_convection` namelist.

    Returns:
        Temperature of *LCL* in *K*.
    """
    def lcl_opt_func(temp_lcl, p_start, temp_start, sphum_start, p_ref):
        # Function to optimize
        r = mixing_ratio_from_sphum(sphum_start)
        theta = potential_temp(temp_start, p_start, p_ref)
        # theta = temp_start * (p_ref / p_start) ** kappa
        value = theta ** (-1 / kappa) * p_ref * r / (epsilon + r)

        return value - saturation_vapor_pressure(temp_lcl) / temp_lcl ** (1 / kappa)
    return optimize.newton(lcl_opt_func, 270, args=(p_start, temp_start, sphum_start, p_ref))

`ref_temp_above_lcl(temp_lcl, p_lcl, p_full)`

Function to replicate the way the reference temperature profile above the LCL is computed in Isca with the Simple Betts-Miller convection scheme.

Parameters:

Name	Type	Description	Default
`temp_lcl`	`float`	Temperature at the Lifting Condensation Level (LCL) in K.	required
`p_lcl`	`float`	Pressure of the LCL in Pa. This will not be one of the pressure levels in the model, and will be larger than any value in `p_full`.	required
`p_full`	`ndarray`	`float [n_p_levels]`. Full model pressure levels in ascending order. `p_full[0]` represents space and `p_full[-1]` is the smallest pressure level in the model that is above the LCL. Units: Pa.	required

Returns:

Type	Description
`ndarray`	`float [n_p_levels]`. Reference temperature in K at each pressure level indicated in `p_full`.

Source code in isca_tools/convection/simple_betts_miller.py

def ref_temp_above_lcl(temp_lcl: float, p_lcl: float, p_full: np.ndarray) -> np.ndarray:
    """
    Function to replicate the way the
    [reference temperature profile above the LCL](https://github.com/ExeClim/Isca/blob/
    7acc5d2c10bfa8f116d9e0f90d535a3067f898cd/src/atmos_param/qe_moist_convection/qe_moist_convection.F90#L621C7-L651)
    is computed in *Isca* with the
    [Simple Betts-Miller](https://jduffield65.github.io/Isca/namelists/convection/qe_moist_convection/)
    convection scheme.

    Args:
        temp_lcl: Temperature at the Lifting Condensation Level (LCL) in *K*.
        p_lcl: Pressure of the LCL in *Pa*. This will not be one of the pressure levels in the model, and will
            be larger than any value in `p_full`.
        p_full: `float [n_p_levels]`.</br>
            Full model pressure levels in ascending order. `p_full[0]` represents space and `p_full[-1]` is the smallest
            pressure level in the model that is above the LCL.</br>Units: *Pa*.

    Returns:
        `float [n_p_levels]`.</br>
            Reference temperature in *K* at each pressure level indicated in `p_full`.
    """
    if np.ediff1d(p_full).min() < 0:
        raise ValueError('p_full needs to be in ascending order.')
    if p_lcl < p_full[-1]:
        raise ValueError(f'p_lcl ({p_lcl}Pa) should be larger than the last value in p_full ({p_full[-1]}Pa).')

    # Add LCL pressure to end of pressure array as it is the largest
    p_full = np.concatenate((p_full, [p_lcl]))
    temp_ref = np.zeros_like(p_full)
    temp_ref[-1] = temp_lcl
    mix_ratio_ref = np.zeros_like(p_full)
    mix_ratio_lcl = mixing_ratio_from_partial_pressure(saturation_vapor_pressure(temp_lcl), p_lcl)
    mix_ratio_ref[-1] = mix_ratio_lcl
    # pressure is in ascending order so iterate from largest pressure (last index) which is LCL, to lowest
    # (first index) which is space
    for k in range(len(p_full) - 2, -1, -1):
        # First get estimate for temperature and mixing ratio half-way between the pressure levels
        # using larger pressure level
        a = kappa * temp_ref[k+1] + L_v/c_p * mix_ratio_ref[k+1]
        b = L_v**2 * mix_ratio_ref[k+1] / (c_p * R_v * temp_ref[k+1]**2)
        dtlnp = a/(1+b)
        temp_half = temp_ref[k+1] + dtlnp * np.log(p_full[k] / p_full[k + 1]) / 2
        mix_ratio_half = mixing_ratio_from_partial_pressure(saturation_vapor_pressure(temp_half),
                                                            (p_full[k] + p_full[k + 1]) / 2)

        # Use halfway values to compute temperature at smaller pressure level
        a = kappa * temp_half + L_v/c_p * mix_ratio_half
        b = L_v**2 * mix_ratio_half / (c_p * R_v * temp_half**2)
        dtlnp = a/(1+b)
        temp_ref[k] = temp_ref[k+1] + dtlnp * np.log(p_full[k] / p_full[k+1])

        # Use temperature at smaller pressure to compute mixing ratio at smaller pressure
        mix_ratio_ref[k] = mixing_ratio_from_partial_pressure(temp_ref[k], p_full[k])
    return temp_ref[:-1]        # don't return LCL value