From one grid point of ERA5 reanalysis pressure-level data, both surface_based_cape_cin() and most_unstable_cape_cin() functions returned a very large negative CAPE. Tried but haven't figured out where a bug might be located in metpy codes.
Here are the example data and code:
import numpy as np
from metpy.units import units
import metpy.calc as mpcalc
pres = np.array([1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1])
temp = np.array([27.8, 26.5, 24.2, 22.1, 20.0, 17.4, 15.2, 13.2, 11.4, 9.6, 7.8, 4.3, 0.5, -3.1, -7.0, -11.5, -16.2, -22.2, -30.3, -39.3, -49.3, -53.7, -57.7, -59.9, -59.7, -60.3, -61.7, -61.3, -57.8, -51.2, -48.2, -37.6, -32.5, -25.9, -16.3, -3.5, -4.9])
tdew = np.array([14.2, 13.2, 12.6, 11.9, 11.2, 10.7, 9.8, 9.1, 7.6, 5.1, 2.0, -2.3, -5.0, -10.6, -19.9, -23.5, -28.2, -32.8, -38.3, -48.7, -58.5, -62.2, -63.1, -68.3, -76.8, -83.7, -85.5, -87.2, -88.7, -91.3, -93.3, -96.6, -98.2, -99.7, -102.1, -103.8, -106.8])
pres = pres * units('hPa')
temp = temp * units('degC')
tdew = tdew * units('degC')
# CAPE/CIN
cape, cin = mpcalc.most_unstable_cape_cin(pres, temp, tdew)
print(f"Most Unstable CAPE: {cape:.1f}, CIN: {cin:.1f}")
# Surface based
cape, cin = mpcalc.surface_based_cape_cin(pres, temp, tdew)
print(f"Surface-based CAPE: {cape:.1f}, CIN: {cin:.1f}")
The output are: Most Unstable CAPE: -109521.2 joule / kilogram, CIN: -65.6 joule / kilogram Surface-based CAPE: -109521.2 joule / kilogram, CIN: -65.6 joule / kilogram
I have tried to manually calculate like the following and obtained a seemingly reasonable CAPE. I also attempted to check some of the metpy source codes. But not sure if there is a bug or not. It seemed that the metpy cape_cin function tried to find an EL and somehow it got Nan and then took the top level (1 mb in the data), which caused the integrated CAPE to be a large negative number. Not sure if find_intersections() had impacted.
# Parcel profile
prof = mpcalc.parcel_profile(pres, temp[0], tdew[0]).to('degC')
print("profile:", "[" + ", ".join(f"{x:.1f}" for x in prof.magnitude) + "]", prof.units)
# Find LFC
lfc_pres, lfc_temp = mpcalc.lfc(pres, temp, tdew, parcel_temperature_profile=prof, which='bottom')
print(f"LFC Pressure: {lfc_pres:.1f}, LFC Temperature: {lfc_temp:.1f}")
# Find EL
el_pres, el_temp = mpcalc.el(pres, temp, tdew, parcel_temperature_profile=prof, which='top')
print(f" EL Pressure: {el_pres:.1f}, EL Temperature: {el_temp:.1f}")
# CAPE calculation between LFC and EL
if lfc_pres is not np.nan and el_pres is not np.nan:
# Find indices for LFC and EL
lfc_idx = np.where(pres >= lfc_pres)[0][-1]
el_idx = np.where(pres <= el_pres)[0][0]
# Integrate CAPE (positive area only)
cape = 0.0
for i in range(lfc_idx, el_idx):
if prof[i] > temp[i]:
delta_ln_p = np.log(pres[i+1].magnitude / pres[i].magnitude)
T_env = temp[i].to('K').magnitude # In Kelvin
T_diff = prof[i].magnitude - temp[i].magnitude
cape += 9.8 * (T_diff / T_env) * delta_ln_p
cape = (-cape * 287 * T_env) * units('J/kg')
else:
cape = 0 * units('J/kg')
# CIN from MetPy
cin = mpcalc.cape_cin(pres, temp, tdew, prof)[1]
print(f"custom CAPE: {cape:.1f}, CIN: {cin:.1f}")
======== output =======
profile: [27.8, 25.6, 23.4, 21.2, 18.9, 16.5, 14.1, 11.7, 10.2, 8.9, 7.6, 4.8, 1.6, -2.0, -6.0, -10.6, -16.1, -22.5, -30.3, -39.7, -50.9, -57.3, -64.4, -72.2, -80.8, -90.5, -101.8, -118.4, -132.6, -151.7, -165.0, -184.4, -193.0, -200.4, -210.2, -217.1, -227.2] degree_Celsius
LFC Pressure: 734.9 hectopascal, LFC Temperature: 6.8 degree_Celsius
EL Pressure: 349.0 hectopascal, EL Temperature: -30.5 degree_Celsius
custom CAPE: 1072.7 joule / kilogram, CIN: -66.3 joule / kilogram
An answer has been found for the causes by @kgoebber and posted at https://github.com/Unidata/MetPy/issues/3790#issuecomment-2807622049.
The cause for both the unreasonably large negative and positive CAPE values were from the mixing-ratio calculation at levels above at least 10 hPa levels, where normal sounding would not be able to reach. By limiting the calculations to levels below 100 hPa (kgoebber) or 50 hPa (my test), the unreasonable values disappear.