most_unstable_cape_cin throwing negative CAPE values

[am-thyst]

I seem to be having problems with the most_unstable_cape_cin function. At the moment, I'm catching the IndexError (#1011 ) in a try-except block and setting these values to zero. But I seem to be getting a high number of negative CAPE values.

This can be shown in the image below, where all the uncoloured areas are negative CAPE:

As far as I can tell, CIN is calculated separately so there should be no negative CAPE. I'm looping over every grid-point in the world, as you can see there are some gridpoints which are being calculated okay.

I chose a grid-point that was bugging as my most negative CAPE value (-657 J/kg). When I plot the profile into a skew-T and tell it to shade CAPE, nothing is shaded. There is no CAPE present:

To replicate my results, here are my profiles and methods:

Td = [293.27972, 289.81912, 287.24918, 281.4828, 269.6113, 270.9546, 273.1311, 274.63028, 265.2233, 194.93207, 187.52663, 183.75009, 181.79565, 180.52945, 176.32547, 171.53108, 166.62222, 160.64876, 154.42905, 152.98645] * units.kelvin

T = [296.38864, 292.4427, 289.31082, 290.3095, 290.9022, 288.90042, 285.52344, 281.09085, 276.7057, 272.3691, 267.90158, 262.1165, 255.22964, 247.16406, 237.98503, 227.96094, 217.04095, 207.87036, 197.53433, 205.37505] * units.kelvin

p = [1009.84, 950, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100, 50] * units.hPa

# CAPE calc: results in -657 J/kg CAPE and -5 J/kg CIN
cape, cin = most_unstable_cape_cin(p, T, Td)

# Skew-T: I usually convert T and Td to degC for this 
fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig)
skew.plot(p, T, 'r', linewidth=2)
skew.plot(p, Td, 'g', linewidth=2)
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
skew.plot(p, parcel_prof, 'k', linewidth=2)
skew.shade_cape(p, t, parcel_prof)
plt.show()

[zbruick]

Thanks for the reproducible example! I just ran this against v0.10.2 and was able to reproduce your error. However, when I run it against master, I get the following output:

print(cape, cin)

(7.977228881942766 <Unit('joule / kilogram')>,
 -675.0915296399907 <Unit('joule / kilogram')>)

This makes a lot more sense, as the large negative number is attributed to CIN, as it should be per the SkewT. I think somewhere in the numerous fixes that will be included in v0.11, this issue was fixed. v0.11 should be out by the end of the month. In the meantime though, if you install MetPy from master, this issue should be resolved. It would be great to have you test this out against your gridded data to make sure!

[am-thyst]

@zbruick with a quick test, it looks fixed. I'll just run this through the entire grid :)

[zbruick]

Ok thanks, keep us posted on your findings!

[am-thyst]

hi @zbruick we reran this across the whole grid and it solved the majority of the problem. There are still some (<0.01%) gridpoints that are throwing this problem. Although it's no longer majorly problematic for us, I thought I should bring it to your attention.

The following profile gives out -27382 J/kg CAPE and -8.918405025406972 J/kg CIN:

Td = [264.5351, 261.13443, 259.0122, 252.30063, 248.58017, 242.66582, 239.64774, 236.17035, 233.02997, 179.90454, 176.28152, 173.54251, 170.59572, 167.68713, 165.8778, 164.73347, 161.91911, 158.15884, 156.13792, 153.60583] # Kelvin

T = [273.09723, 268.40173, 263.56207, 260.257, 256.63538, 252.91345, 248.48195, 243.75995, 239.11597, 235.1402, 231.77188, 227.87416, 223.00409, 216.6683, 210.2053, 206.43655, 207.02776, 206.75017, 203.48286, 197.07687] # Kelvin

p = [1017.16, 950, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100, 50] # hPa

and the following skew-T:

[zbruick]

Haha ouch. That's not great. Thanks for providing this profile. We'll look into it and see what's going on.

[am-thyst]

@zbruick I was wrong saying it was only <0.01% points throwing negative values, I think I interpreted my colleague's message wrong. The upgrade to master has only fixed a few points. The mid-latitudes and high latitudes are still bugging quite severely.

Here's a before and after:
Before:

After:

[zbruick]

So I think what's happening here (and resultantly what the change in #1172 fixes) is that the model you're using has a pretty low vertical resolution. As a result, it can easily happen that there is positive area below the LCL, meaning that the LFC and EL both occur below the LCL. Our current code did not catch this exception before, meaning that the LFC was bumped to the LCL and the EL was never found as it didn't occur above the LCL/LFC. So CAPE was calculated through the full profile, which was all negative area. The change in #1172 will make the LFC and EL NaN values, which will prevent CAPE from being calculated (it'll return NaN), which is the appropriate answer, at least for the profile you provided above.

[zbruick]

If you want, you can clone my repo and pull the branch used in #1172 to test to see if that does work for you.

[am-thyst]

@zbruick fixed! Thanks