-
Notifications
You must be signed in to change notification settings - Fork 415
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add downdraft CAPE calculation #823
Conversation
depth=pressure[-1] - top_pressure) | ||
|
||
|
||
# Calculate the difference in profile and environmental temperature to integrate |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
E303 too many blank lines (2)
Ok working on this some more there are some implementation details to think about. The kwargs |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It looks liketop_pressure
and bottom
are actually used when parcel
is passed in.
metpy/calc/thermo.py
Outdated
dewpoint = dewpoint[sort_inds] | ||
|
||
# Trim data to only top height and below as well as bottom height and above | ||
pressure, temperature, dewpoint = get_layer(pressure, temperature, dewpoint, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This call of get_layer
only differs from the one below with the passing of bottom
. Are they both needed?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That was a try at something that won't survive 😄 I think I'm going to refactor some here - we don't want to trim before we interpolate a parcel.... hmmm.. more to come.
metpy/calc/thermo.py
Outdated
|
||
# Trim data to parcel starting point and below, interpolating the starting point | ||
pressure, temperature, dewpoint = get_layer(pressure, temperature, dewpoint, | ||
depth=np.nanmax(pressure) * pressure.units - parcel_starting_pressure, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
E501 line too long (122 > 95 characters)
heights=heights) | ||
|
||
|
||
# The user did not give us a parcel, so we'll calculate a sensible default. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
E303 too many blank lines (2)
metpy/calc/thermo.py
Outdated
|
||
# Trim data to only top height and below as well as bottom height and above | ||
pressure, temperature, dewpoint = get_layer(pressure, temperature, dewpoint, | ||
depth=np.nanmax(pressure) * pressure.units - top_pressure, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
E501 line too long (110 > 95 characters)
parcel=custom_parcel) | ||
dcape_truth = 101.56717405359117 * units.joule / units.kilogram | ||
dcape_pressure_truth = np.array([973, 943.5, 925, 910.6, 878.4, 865, 850, 848, 847.2, | ||
816.9, 793, 787.8, 759.6, 759, 732.2, 700, 670]) * units.hPa |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
E501 line too long (97 > 95 characters)
metpy/calc/thermo.py
Outdated
|
||
# Trim data to parcel starting point and below, interpolating the starting point | ||
pressure, temperature = get_layer(pressure, temperature, | ||
depth= profile_depth, heights=heights) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
E251 unexpected spaces around keyword / parameter equals
metpy/calc/thermo.py
Outdated
# Trim data to be only above the bottom, interpolating if necessary. This must be done | ||
# in separate step as there is no way to specify the top of a layer, just depth. | ||
pressure, temperature = get_layer(pressure, temperature, bottom=bottom, | ||
depth= profile_depth, heights=heights) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
E251 unexpected spaces around keyword / parameter equals
So the RAOB program also calculates DCAPE and can do so in three different ways. It appears to do it via a lowest depth (e.g., 6 km) through either a density weighted average (most accurate), average wet-bulb, or coldest wet bulb (least accurate) [accuracy info from RAOB program]. Attached are two images using the coldest wet bulb temperature over the lowest 6 km from RAOB for the two times in the above examples from OUN for 17 April 2018 (ignore the month - its the same data). This gives different values from the NSHARP values given in the above figures - likely a different method. |
Speaking of the NWS method (is this the NSHARP way?), I found the following PDF that talks about the DCAPE process (see slide number 6 DCAPE_Web.pdf). To summarize the process:
|
Are there any updates to the status of this, or what needs to be done before this can get merged? This would be very helpful for a couple ongoing projects in the research group I'm in. |
@jthielen So we just need some help validating that what we're doing is actually correct. A scholarly source to cite for the method plus some good values to double check our work are what's needed. |
I'm working in a proyect and I need the calculation of the Downdraft CAPE (DCAPE) and I'm testing this function and the other mentioned here. My surprise is that I found diferent results if I compared whith the data of the sounding. I choose this sounding (https://www.spc.noaa.gov/exper/soundings/22081900_OBS/) and I realized the calculus whith the function downdraft_cape described in this post and the function created by jillianndufort and described in https://github.com/jillianndufort/Calculation-for-DCAPE I found different results for the same data and differents of the result DCAPE in the sounding. 1.- This function downdraft_cape: I don't understand what is happening because I think that the functions are OK in the sense of the programming. I attached the data file so you can do the tests. One last comment, in this downdraft cape function, I found a problem in the dcape.to('J/kg') conversion giving me an error: DimensionalityError: Cannot convert from 'delta_degree_Celsius * joule / kilogram' ([length] ** 2 * [temperature] / [time] ** 2) to 'joule / kilogram' ([length] ** 2 / [time] ** 2) Ignoring this, the result is the commented. Could you help me? |
Thanks for kicking the tires on this calculation. Given that unit error, the implementation of the calculation is definitely not correct. Whether that's causing the difference in values, though, is unclear. Unfortunately, I don't have any cycles at the moment to dig into this. |
I also think that they are two independent problems: the problem with the units and the calculus of the DCAPE. I'll try to check this these days, before my holidays. I think that it doesn't seem easy to solve. |
Superseded by #3120. |
Closes #728 by adding a downdraft CAPE calculation. As it stands there are still some things to discuss/clear up:
The main difficulty currently is that this does not appear to compare well with the SPC DCAPE values. Looking at their results I'm not sure there isn't a bug on one or both ends. Take the two soundings below. To me one should have positive DCAPE, the other negative as in one the parcel path is to the left of the environmental temperature, and it is to the right in the other. Also, to get values anywhere near as large as these I have to integrate the entire lower 400 hPa, not the surface up to the minimum theta-e in the lowest 400 hPa. Any comments welcome as we need to nail this down before testing and merging.