Lifting Parcels with SHARPlib
In this article we’re going to take a look at how to use SHARPlib for a fairly common use case- computing parcel-based moist adiabatic ascent. Many parts of SHARPlib are familiar in nature to SHARPpy and NSHARP, but this is an area where SHARPlib can differ significantly. This was a conscious choice in order to generalize the process of computing parcel-based ascent so that different methods of computing said ascent can be used, tested, and configured for specific applications. This example should also be helpful in providing an example of the sort of preprocessing and formatting the library expects from the programmer/user before being called. This article is primarily going to focus on using the Python wrapped version of the library - if you’re interested in a C++ demo, checkout the GitHub reporistory folder located at: SHARPlib/examples/C++.
Step 0: Imports, I/O, and Preprocessing
SHARPlib has a few basic requirements when working with it:
- 1) You handle your own reading of data (text, grib2, NetCDF)
- 2) Profile data must be represented as 1D arrays
- 3) All units are assumed to be base MKS (meters, kelvin, seconds, unitless mixing ratio, etc)
Where SHARPpy hid some of these extra steps behind the scenes to leverage creative reuse of computations, SHARPlib will require the programmer to be explicit about what they want. The library tries its best to abide by a “don’t pay for what you don’t use” philosophy. However, these operations aren’t too difficult: an example of reading a sounding file, converting units, and computing additional arrays of profile data are shown below. Of particular importance are the creation of virtual temperature and mixing ratio arrays, which will be needed for parcel ascent.
from nwsspc.sharp.calc import layer
from nwsspc.sharp.calc import parcel
from nwsspc.sharp.calc import interp
from nwsspc.sharp.calc import thermo
from nwsspc.sharp.calc import constants
import pandas as pd
filename = "https://raw.githubusercontent.com/keltonhalbert/SHARPlib/main/data/test_snds/20160524_2302_EF3_37.57_-100.13_108_613967.snd"
col_names = ["pres", "hght", "tmpc", "dwpc", "wdir", "wspd", "omeg"]
snd_data = pd.read_csv(filename, skiprows=7, comment="%", names=col_names)
## preprocess additional fields
snd_data["pres"] = snd_data["pres"] * constants.HPA_TO_PA
snd_data["tmpk"] = snd_data["tmpc"] + constants.ZEROCNK
snd_data["dwpk"] = snd_data["dwpc"] + constants.ZEROCNK
snd_data["mixr"] = thermo.mixratio(snd_data["pres"], snd_data["dwpk"])
snd_data["vtmp"] = thermo.virtual_temperature(snd_data["tmpk"], snd_data["mixr"])
snd_data["thta"] = thermo.theta(snd_data["pres"], snd_data["tmpk"])You should be able to copy this code verbatim into a script or Jupyter notebook, and it will read the remote file in the GitHub repository, parse it, and compute the requested fields. This particular sounding file comes from the SPC Mesoanalysis, but the process is similar for both observed and gridded sounding data. I’ll discuss it in a future blog post, but if you do want to attempt to work with a gridded dataset, apply this preprocessing step to each vertical profile on the grid and construct a 1D flat array of profiles of size NX*NY! For now, however, we’ll stick to a single profile. Now that we’ve read and preprocessed our data, we can move on to the more interesting stuff.
Step 1: Define Your Parcel and Moist Adiabat
To compute moist adiabatic ascent, we need to configure two things:
- 1) Our air parcel’s starting pressure, temperature, and dewpoint temperature
- 2) Which solver to use for moist ascent, and its associated configuration options
For simplicity in this first example, we will use the Wobus pseudoadiabat used by NSHARP and SHARPpy in operational environments. It’s known to have errors in parcel temperatures of up to 1K when compared to more accurate solvers, but is very fast. It also doesn’t require any additional configuration to work with, meaning we can just use it right out of the box! This technically uses a polynomial fit to the Wobus function, which is defined by Table 78 of PP.319-322 of the Smithsonian Meteorological Tables by Roland List (6th Revised Edition). More information about the errors with the Wobus function can be found in Davies-Jones 2008.
The simplest parcel we can define is the surface based parcel, where we set the parcel starting attributes to the values of the environment near the surface (in this case, index 0). Later we will discuss how to define something more complex, such as a mixed-layer parcel.
## Instantiate a parcel lifter using
## the Wobus liquid-only pseudoadiabat.
## This has been the operational default for
## NSHARP/SHARPpy due to its speed, but
## is known to be somewhat inaccurate.
lft_wobf = parcel.lifter_wobus()
## define a surface based parcel by instantiating
## a parcel object, and setting it's initial
## pressure (Pa), temperature (K), and dewpoint (K).
sb_pcl = parcel.Parcel()
sb_pcl.pres = snd_data["pres"].iloc[0]
sb_pcl.tmpk = snd_data["tmpk"].iloc[0]
sb_pcl.dwpk = snd_data["dwpk"].iloc[0]Step 2: Compute and Integrate buoyancy
Now that we’ve defined our parcel starting attributes, and we have chosen a solver for moist adiabatic ascent, we can compute buoyancy. Once we have buoyancy, we can integrate it to get CAPE and CINH! All buoyancy calculations in SHARPlib use virtual temperature (also known as the virtual temperature correction), since the presence of water vapor can have a significant enough impact in small CAPE environments. If you’d like to learn more about the impacts of neglecting the virtual temperature correction, check out Doswell and Rasmussen 1994. Calling lift_parcel will also set the pressure level of the LCL in sb_pcl.lcl_pressure.
Finally, we can integrate the buoyancy to get CAPE and CINH. The cape_cinh routine will first call the find_lfc_el algorithm in order to determine the integration bounds. For SPC use cases, we’ve opted for the LFC-EL combination that results in the maximum buoyancy for the profile, since we’re mostly concerned with convection. However, it is possible to manually define integration intervals, which will be discussed further below.
## lift_parcel sets the sb_pcl.lcl_pressure
## attribute for us, in addition to returning
## an array of buoyancy
buoyancy = parcel.lift_parcel(
lft_wobf,
snd_data["pres"],
snd_data["vtmp"],
sb_pcl
)
## cape_cinh calls find_lfc_el for us using
## an algorithm that searches for the LFC-EL
## layer with the maximum buoyancy.
sb_cape, sb_cinh = parcel.cape_cinh(
snd_data["pres"],
snd_data["hght"],
buoyancy,
sb_pcl
)For the simplest example, that’s all there is to it! While other details may differ depending on the type of parcel you are lifting, or the type of moist ascent solver you use, the process is largely the same. It should be noted, however, that this API is still a work in progress and may be subject to change, especially if I get any particularly insightful feedback from users on this process. In the absence of that, I only expect the function names and call signatures to change slightly, but not dramatically. In SHARPpy and NSHARP, a singular routine was responsible for most of these operations. By splitting them out into their component parts, it allows for more general code reuse and flexibility of methods. I am personally biased towards this approach because to me, it comes across as a more process or ingredients-based approach that flows very naturally if you understand the mechanics of parcel lifting. If you don’t fully understand the mechanics of parcel lifting, then the hope is that constructing the API in this way makes it clear what is going on in the process.
These are just the very basics of what you can do, however. Below are some other related tasks and operations relevant to parcel lifting that may be of interest.
Q: Calling lift_parcel will return buoyancy, but what if I want to display the parcel virtual temperature trace?
There is currently no explicit function that returns only the virtual temperature trace of parcel ascent. At the time of writing the API, my impression was that the vast majority of the time, the user will want buoyancy rather than the raw temperature trace. However, in the process of writing this tutorial, I wonder if that was perhaps an oversight. I would love feedback on what the desired behavior should be, if you have any! Still, it is possible to get a parcel virtual temperature trace, and it also provides a nice means of explaining what happens under the hood when you call lift_parcel - because these operations are effectively identical, sans the computation of buoyancy.
This code block assumes you have already run Step 0 and Step 1. We first compute the LCL pressure and temperature, and assign the LCL pressure to the surface based parcel. Then we define layers for dry and saturated ascent, keep virtual potential temperature constant below the LCL, and call our moist adiabat function to get the virtual temperature of the parcel as it ascends the profile.
## Manually compute the LCL and set the pressure level of the
## LCL to the surface parcel.
lcl_p, lcl_t = thermo.drylift(sb_pcl.pres, sb_pcl.tmpk, sb_pcl.dwpk)
sb_pcl.lcl_pressure = lcl_p
## We have two layers of ascent: one for dry ascent, from the parcel
## starting pressure until the LCL, and another from saturated ascent,
## defined from the lcl pressure until the end of the profile.
dry_lyr = layer.PressureLayer(sb_pcl.pres, sb_pcl.lcl_pressure)
sat_lyr = layer.PressureLayer(sb_pcl.lcl_pressure, snd_data["pres"].iloc[-1])
## We want the array indices that correspond to these layers
dry_idx = layer.get_layer_index(dry_lyr, snd_data["pres"])
sat_idx = layer.get_layer_index(sat_lyr, snd_data["pres"])
## Parcel virtual potential temperature is conserved
## from the initial lift until the LCL is reached.
thetav_lcl = thermo.theta(sb_pcl.lcl_pressure, lcl_t, constants.THETA_REF_PRESSURE)
thetav_lcl = thermo.virtual_temperature(thetav_lcl, thermo.mixratio(sb_pcl.pres, sb_pcl.dwpk))
## Allocate space for our output arrays.
tmpk_arr = np.zeros(snd_data["pres"].shape, dtype="float32")
pres_arr = np.zeros(snd_data["pres"].shape, dtype="float32")
## dry ascent - we go to ktop+1 because we want to include
## the value of ktop in the loop.
for idx in range(dry_idx.kbot, dry_idx.ktop+1):
tmpk_arr[idx] = thermo.theta(constants.THETA_REF_PRESSURE, thetav_lcl, snd_data["pres"].iloc[idx])
## saturated ascent- we go to ktop+1 because we want to include
## the value of ktop in the loop.
for idx in range(sat_idx.kbot, sat_idx.ktop+1):
tmpk_arr[idx] = lft_wobf(sb_pcl.lcl_pressure, lcl_t, snd_data["pres"].iloc[idx])This is likely to come across as overly involved and verbose - but as mentioned above, if there’s enough interest, this could be standardized into the API rather easily. Part of the goal in rewriting this library was to also provide a set of tools with which to develop new algorithms and techniques. So, while this is meant to showcase how to get a parcel trace, it also showcases some of the individual underlying components and how they can be composed into interesting operations.
Q: What if I want to compute a mixed-layer parcel?
Since we manually set parcel starting attributes in Step 1, it’s fairly straightforward to define something like a mixed-layer parcel (and even use your own mixing depth). Execute Step 0 and then replace Step 1 with the following code…
## Same as before in Step 1
lft_wobf = parcel.lifter_wobus()
mix_depth = 10000.0 ## 100 hPa depth
sfc_pres = snd_data["pres"].iloc[0]
mix_lyr = layer.PressureLayer(sfc_pres, sfc_pres - mix_depth)
## compute the mean theta and mixing ratio over the layer
mean_mixr = layer.layer_mean(mix_lyr, snd_data["pres"], snd_data["mixr"])
mean_thta = layer.layer_mean(mix_lyr, snd_data["pres"], snd_data["thta"])
## define a surface based parcel by instantiating
## a parcel object, and setting it's initial
## pressure (Pa), temperature (K), and dewpoint (K).
ml_pcl = parcel.Parcel()
ml_pcl.pres = snd_data["pres"].iloc[0]
## invert theta for the starting temperature
ml_pcl.tmpk = thermo.theta(constants.THETA_REF_PRESSURE, mean_thta, sfc_pres)
## get the dewpoint from mixing ratio
ml_pcl.dwpk = thermo.temperature_at_mixratio(mean_mixr, sfc_pres)Then, you may proceed with Step 2 in order to get buoyancy, CAPE, and CINH.
Q: What if I only want to integrate buoyancy over some arbitrary layer, i.e. 0-3km CAPE?
One of the many code patterns observed when looking at the operational code at SPC was integrating buoyancy over arbitrarily defined layers, i.e. 0-3km CAPE, 0C to -20C CAPE, and others. A revised Step 2 to accomplish this would look like the following code, presuming Step 0 and Step 1 have already been executed.
## lift_parcel sets the sb_pcl.lcl_pressure
## attribute for us, in addition to returning
## an array of buoyancy
buoyancy = parcel.lift_parcel(
lft_wobf,
snd_data["pres"],
snd_data["vtmp"],
sb_pcl
)
## Find the LFC-EL combination with maximum buoyancy, and set
## the values of those pressure levels in sb_pcl.lfc_pressure
## and sb_pcl.eql_pressure.
parcel.find_lfc_el(sb_pcl, snd_data["pres"], snd_data["hght"], buoyancy)
## Define the layer over which to accumulate CINH.
lpl_lfc_plyr = layer.PressureLayer(sb_pcl.pres, sb_pcl.lfc_pressure)
## integration of buoyancy requires height coordinates for proper units,
## so convert the pressure layer to a height layer
lpl_lfc_hlyr = layer.pressure_layer_to_height(lpl_lfc_plyr, snd_data["pres"], snd_data["hght"])
## Create the 0-3km HeightLayer for integration.
sfc_hght = snd_data["hght"].iloc[0]
hlyr_0_3km_agl = layer.HeightLayer(sfc_hght, sfc_hght + 3000.0)
## Integrate the buoyancy using the trapezoid method. The -1 argument signals the intent to
## only integrate negative area, while the +1 argument signals the intent to integrate
## positive area. An argument of 0 would integrate both positive and negative values.
sb_pcl.cinh = layer.integrate_layer_trapz(lpl_lfc_hlyr, buoyancy, snd_data["hght"] -1);
sb_pcl.cape = layer.integrate_layer_trapz(hlyr_0_3km_agl, buoyancy, snd_data["hght"], 1);Using this approach, it becomes pretty straightforward to define any arbitrary interval over which to integrate buoyancy.
Q: How do I configure a more accurate moist adiabat?
The Wobus pseudoadiabat is nice and fast, but it only accounts for liquid hydrometeors and does not include any latent heat release from ice hydrometeors. As of the writing of this article, the only other supported solver for moist ascent is lifter_cm1, derived from George Bryan’s CM1 CAPE routine and the Bryan and Fritch 2004 paper on ice-liquid water potential temperature. The current implementation of this algorithm is a little less performant than I know it can be, so it may start to bottleneck on gridded datasets until I can make the appropriate revisions. Still, this should give you a comprehensive idea on how to configure it. Execute Step 0, and then supplement Step 1 with the code below…
lft_cm1 = parcel.lifter_cm1()
## For this example, we will use a pseudoadiabat
## (does not retain liquid during ascent) that includes
## ice processes
lft_cm1.ma_type = thermo.adiabat_pseudo_ice
## The solver integrates in defined pressure increments,
## and for accuracy, we will set the interval to 1 hPa.
## You may want to experiment to find an optimal tradeoff
## between speed and accuracy.
lft_cm1.pressure_incr = 100.0
## The convergence criteria for the iterative solver. In this case,
## the output must be within 0.001 Kelvin of the solution to count
## as converged.
lft_cm1.converge = 0.001;Once you’ve configured the parcel lifter, you may proceed with Step 2 as normal. The CM1 parcel lifter can be configured for 4 types of saturated ascent:
- Pseudoadiabatic ascent, liquid only
- Pseudoadiabatic ascent, liquid + ice
- Adiabatic ascent, liquid only
- Adiabatic ascent, liquid + ice
The adiabatic ascent is analogous to what some would call “reversible CAPE”. However, I’ve been corrected in this a number of times that it is technically not reversible since entropy still increases, but the type of ascent is still adiabatic in nature.
Wrapping it all up…
Hopefully this has been helpful in showcasing how the logic and flow of working with SHARPlib was intended to work. However, what seems intuitive, clear, and readable to me may not be clear and readable to everyone. The advantage of releasing SHARPlib in its current “early access” period is that it means I get to query you all on usage, points of friction, and points of confusion! So, if you have any feedback (good or bad), feel free to shoot me a message via email or on your platform of choice. Ideally, I’d also like to convert these blog posts into a Jupyter notebook or series of notebooks - perhaps using Project Pythia - as a means of providing interactive documentation and examples. For now, however, while things are still fluid and in development, it’s convenient to be able to dump things here on the blog!
Happy parcel lifting,
Kelton.