Upgrading parcel lifting routines

[wblumberg]

So, this issue is mainly a TODO list for myself and to announce some of the future plans for upgrading SHARPpy's parcel lifting routines. One of the key things that has been said by users is that the parcel lifting routines tend to be slow. I've been working for various ways to get around this:

1.) Upgrade the method of calculating pseudoadiabats (or saturated parcels in general) by using the Davies Jones (2008) method. I've been experimenting with this and have incorporated it into the parcelx() and cape() routines used by SHARPpy to do parcel lifting type calculations. While this method has not sped up the routines significantly when using the typical standard and significant level calculations (~200 points), but when testing full resolution soundings, the speed up is roughly 2x.

2.) This upgrade has facilitated revisiting of the logic of the effective_inflow_layer() and convective_temp() functions, which are computationally expensive to calculate. One of the ideas I've tossed around is taking advantage of the fact that when lifting parcels to determine if CAPE or CIN meets a certain threshold, you do not have to iterate throughout the entire profile. The loop can be broken once the threshold is met, which will speed up the overall routine. For the effective_inflow_layer(), the checks for the CAPE and CIN meeting thresholds can be switched on and off. The AND and the OR operators can also be switched to break the loop during the lifting process.

3.) The cape_rdj() and convective_temp calculations could also have pseudoadiabats calculated en mass prior to the actual call to the parcel lifting (the RDJ method enables that function).

[wblumberg]

I've successfully implemented logic for points 1 and 2. The effective_inflow_layer calculation is definitely faster...the question is by how much. I can't really think of anyway to improve the convective_temp() function yet using the changes I've made. Maybe there's a better search algorithm for to search for a parcel that has CINH = 0 J/kg than the current method. 3 may wait, since it may be a little more intensive.

TODO: Fix thermo.wetbulb() to support both Wobus and RDJ-2008 methods of calculating wet bulb temperature (and therefore the wet bulb temperature profile. TODO: Modify the dcape() routine to provide the option of using the Wobus vs. RDJ 2008 methods. TODO: Modify the wet bulb precipitation type algorithm [watch_type.posneg_wetbulb()] to support using the Wobus vs RDJ 2008 methods.

Do timing and value comparison tests using the SARS soundings. Compare:

CAPE/CIN/LCL/LFC/EL Effective inflow layer top and bottom DCAPE wet bulb and theta-e profiles.

[tennessean]

Recently, I found another method for calculating saturated parcels that uses polynomials instead of iterations:

Moisseeva, Nadya and Roland Stull. (2017) Technical note: A noniterative approach to modelling moist thermodynamics. Atmospheric Chemistry and Physics 17:24, 15037-15043.

(The article and its supplement are both freely available for download.)

The article describes the development of a polynomial "best fit" for constructing moist adiabats, going into some detail about how it's constructed. The main calculator (and its governing equations) is provided in an Excel document in the supplement. The equations themselves are fairly complicated and may need some modification from how they're structured in the Excel format, but at the same time the article notes that the calculations themselves are actually less computationally expensive than the iterative methods, and by several orders of magnitude to boot!

In addition, the accuracy of the equations is noted to be to within perhaps a few tenths of a degree Celsius from the graphical method, but they also note that the various types of thermodynamic diagrams (emagram, tephigram, Stuve, Skew-T) have similar slight differences regarding the coordinate positions of their respective moist adiabats.

All in all, I at least think that it should be worth at least taking a look.

Please let me know what you think about this suggestion.

Thank you, and have a good day.

[wblumberg]

@tennessean thank you for the heads up. I have seen this paper before, and I am encouraged by the speed of this method. However, because a lot of the features of SHARPpy (SARS, the Sig Tor/EF climatology) are calibrated using the Wobus function inherent in SHARP, replacing the logic behind calculating the pseudoadiabats may reduce the usefulness of these features. In addition, a motivation to develop SHARPpy was because the variations in the parcel lifting code between different software packages cause noticeable differences in the CAPE and CIN values that were computed. We've seen this when comparing a variety of radiosonde analysis packages, and I've even seen this in my tests using the RDJ parcel lifting code mentioned above. Some of my own work studying how sensitive the methodologies behind the calculation of CAPE and CIN are make me extremely hesitant to make even subtle changes to the code. Therefore, I'm very cautious about replacing key features in this software without serious consideration about how to inform users about its implications as it's antithetical to the reproduction crisis we're trying to help solve. I've noticed that journal articles don't even reference the version number of the SHARPpy code they use, which will impact the reproducibility of science (and there are small differences).

Your recommendation makes me think I'm going to go ahead and code this up to see what it does to my comparisons. If you have any recommendations for how to go about informing users about these changes (or other aspects of this problem), I'm all ears.

[tennessean]

One idea that I came up with was to have the option of a drop-down menu or somesuch that would allow one to select among the various types of moist-adiabat calculation methods (Wobus, Davies-Jones, Moisseeva/Stull, and perhaps even others). This could allow for relatively easy comparison of data, and the selected formula could be subsequently mentioned in articles making use of SHARPpy (e.g., "The moist-adiabat formula that we chose to use was the X formula, for Y reasons.").

Of course, I'm admittedly not much of a programmer, and something like this is beyond my limited skillset. I'm not even sure if it's actually possible. Nonetheless, I would at least hope that some kind of customization of this kind would be workable.