WK82 analytical sounding

[aschueth]

I'm not sure where this would exactly fit in the ecosystem, but I think it might be useful for us cloud modelers to have analytic soundings and wind profiles in metpy. It might be useful for the sounding examples to have this built in, so no text file or downloading sounding data is required. Furthermore, having the code built in could be helpful to have it in case people (like me) want to manipulate it for experiments. This might be a low use-case situation, but it still could be helpful.

Here is a python implementation in another package: https://github.com/cwebster2/pyMeteo/blob/master/pymeteo/cm1/soundings/WK82.py

[dopplershift]

In general, I think functionality like this is completely appropriate for MetPy. I'd have to defer to others' opinions on how generally useful a specific example is. In this case, though, the original paper looks really highly cited. So I'm ok with adding this one.

So then the question becomes: where in MetPy should this live? Maybe just a new module in calc?

Would you be interested in submitting a function that returns the appropriate profile as an xarray Dataset (with units encoded as attribute metadata)?

[aschueth]

I was wondering the same, whether it should live on it's own in calc, or maybe be under the soundings umbrella? If you are feeling really ambitious, it could also be added into some testing, since it should always give constant indices? Here is the function I came up with. I tested it out by plotting the skew-t and it gives expected values. I hope the doc string is thorough enough and in the format you need? I also included the required imports in case they were needed. Is this kind of what you were hoping for? Any input on anything?

import xarray as xr
import metpy.calc as mpcalc
import metpy.constants as mpconst
from metpy.units import units
import numpy as np

def WK82(altitude=np.arange(18000)*units.m, potential_temperature_sfc=300.*units.degK, mixing_ratio_pbl = 0.014*units('kg/kg'), altitude_tropopause = 12000.*units.m, temperature_tropopause=213.*units.degK, potential_temperature_tropopause=343.*units.degK, pressure_sfc=1000*units.hPa):

    r"""Calculate the Weisman and Klemp analytical thermodynamic profile used for idealized cloud models.

    This calculation has the default quantities that can be changes as the arguments in the function. 
    The implementation uses the formula outlined in [WeismanKlemp1982] pg.506.

    Parameters
    ----------
    altitude : `pint.Quantity`
        Heights AGL
    potential_temperature_sfc : `pint.Quantity`
        Potential Temperature at the surface
    mixing_ratio_pbl : `pint.Quantity`
        Constant mixing ratio in the planetary boundary layer
    altitude_tropopause : `pint.Quantity`
        Altitude of tropopause
    temperature_tropopause : `pint.Quantity`
        Temperature at tropopause
    potential_temperature_tropopause : `pint.Quantity`
        Potential temperature at tropopause
    pressure_sfc : `pint.Quantity`
        Pressure at the surface

    Returns
    -------
    Xarray Dataset, attributes contain units

    altitude :
        Heights AGL
    pressure :
        Pressure values at every height level
    potential_temperature : 
        Potential temperature values at every height level
    temperature : 
        Temperature values at every height level
    dewpoint :
        Dewpoint temperature values at every height level
    relative_humidity :
        Relative humidity values at every height level
    mixing_ratio :
        Mixing ratio values at every height level

    """

    potential_temperature = potential_temperature_sfc + (potential_temperature_tropopause - potential_temperature_sfc)*((altitude/altitude_tropopause)**(1.25))
    rh = 1.0-0.75*((altitude/altitude_tropopause)**1.25)

    potential_temperature[int(altitude_tropopause.m):] = potential_temperature_tropopause * np.exp((mpconst.g/(mpconst.Cp_d*temperature_tropopause))*(altitude[int(altitude_tropopause.m):]-altitude_tropopause))
    rh[int(altitude_tropopause.m):] = rh[int(altitude_tropopause.m)]

    pressure = mpcalc.add_height_to_pressure(pressure_sfc, altitude)
    temperature = mpcalc.temperature_from_potential_temperature(pressure,potential_temperature)
    qv = mpcalc.mixing_ratio_from_relative_humidity(rh,temperature,pressure)

    #sets PBL mixing ratio constant up until it is near equal to previous profile
    qv[:(np.abs(qv.m - mixing_ratio_pbl.m)).argmin()]=mixing_ratio_pbl 

    rh=mpcalc.relative_humidity_from_mixing_ratio(qv,temperature,pressure)
    dewpoint=mpcalc.dewpoint_from_relative_humidity(temperature,rh)

    dims = ['nk']
    coords = {
        'altitude': (dims,altitude.m, {'units':altitude.units})
    }

    data_vars = {
        'pressure': (dims, pressure.m, {'units':pressure.units}),
        'potential_temperature': (dims, potential_temperature.m, {'units':potential_temperature.units}),
        'temperature': (dims, temperature.m, {'units':temperature.units}),
        'dewpoint': (dims, dewpoint.m, {'units':dewpoint.units}),
        'relative_humidity': (dims, rh.m, {'units':rh.units}),
        'mixing_ratio': (dims, qv.m, {'units':qv.units})
    }
    ds = xr.Dataset(data_vars, coords)

    return ds

[dopplershift]

Looks great!

I think this belongs in calc along with the sounding indices--then we could document it alongside the rest of the sounding calculations. Docstring looks good. There are a few long lines and missing spaces that will need to be cleaned up to make the bots happy. Would love to see an automated test added in test_indices.py that would check whatever values you checked to validate.

Interested in starting a PR with this and we can finish any edits over there?