Skip to content
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

Closed
wants to merge 9 commits into from
Closed

Conversation

jrleeman
Copy link
Contributor

@jrleeman jrleeman commented Apr 24, 2018

Closes #728 by adding a downdraft CAPE calculation. As it stands there are still some things to discuss/clear up:

  • How should we let folks specify custom bounds for the integration? A custom parcel?
  • Need to add tests

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.
oun-2
oun-3

@jrleeman jrleeman added the Area: Calc Pertains to calculations label Apr 24, 2018
@jrleeman jrleeman added this to the 0.8 milestone Apr 24, 2018
depth=pressure[-1] - top_pressure)


# Calculate the difference in profile and environmental temperature to integrate

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)

@dopplershift dopplershift added the Type: Feature New functionality label May 1, 2018
@jrleeman
Copy link
Contributor Author

jrleeman commented May 1, 2018

Ok working on this some more there are some implementation details to think about. The kwargs bottom and top make sense when using the default parcel. But when using a parcel kwarg that is currently a tuple of pressure and temperature at which we want to start a parcel. It would then go down the moist adiabat to the surface. Then the bottom and top_pressure arguments are ambiguous as to their utilization. Maybe we only take the profile as a complete profile that the user wants? Thoughts @dopplershift?

Copy link
Member

@dopplershift dopplershift left a 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.

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,
Copy link
Member

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?

Copy link
Contributor Author

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.


# 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,

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.

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)


# 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,

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

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)


# Trim data to parcel starting point and below, interpolating the starting point
pressure, temperature = get_layer(pressure, temperature,
depth= profile_depth, heights=heights)

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

# 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)

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

@jrleeman jrleeman modified the milestones: 0.8, 0.9 May 4, 2018
@kgoebber
Copy link
Collaborator

kgoebber commented Jul 9, 2018

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.

ii_jjermq6w1_16480e600c7e473c
ii_jjermq6o0_16480e600c7e473c

@kgoebber
Copy link
Collaborator

kgoebber commented Jul 9, 2018

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:

  1. From 700 mb (or anywhere 500-700 mb), determine the LCL.
  2. Follow moist adiabat to surface from LCL found in 1.
  3. DCAPE is the negative area (positive value) between the downdraft moist adiabat and the temperature profile.

@dopplershift dopplershift modified the milestones: 0.9, 0.10 Jul 26, 2018
@jrleeman jrleeman removed this from the 0.10 milestone Dec 12, 2018
@jthielen
Copy link
Collaborator

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.

@dopplershift
Copy link
Member

@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.

Base automatically changed from master to main February 22, 2021 22:39
@mtesouro
Copy link

mtesouro commented Aug 19, 2022

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:
dcape = 117.76574757857342 delta_degree_Celsius * joule / kilogram
2.- The jillianndofort's function:
{'DD_method1': {'metpy': 41165.230535591356,
'Magnus': 0.0,
'Wexler': 0.0,
'Buck': 0.0},
'DD_method2': {'metpy': 473702.59788943676,
'Magnus': 2.241814634518002,
'Wexler': 2.2423655314996664,
'Buck': 2.242472299362524},
'DD_method3': {'metpy': -71.44700746827606,
'Magnus': 26.86878193998087,
'Wexler': nan,
'Buck': 29.057112541161665}}
3.- Sounding:
DCAPE = 983 j/KG

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.
I have the same problem when I check this function whith the test_dcape_defaults. Whithout to('J/kg') it works fine.

Could you help me?
MFL
datos.txt

@dopplershift
Copy link
Member

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.

@mtesouro
Copy link

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.
If someone can have a look, will be welcome.
Thanks a lot.

@dopplershift dopplershift mentioned this pull request Jul 28, 2023
2 tasks
@dopplershift
Copy link
Member

Superseded by #3120.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Area: Calc Pertains to calculations Type: Feature New functionality
Projects
None yet
Development

Successfully merging this pull request may close these issues.

DCAPE for MetPy
6 participants