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 Galvez-Davison Index calculation #3121

Merged
merged 10 commits into from
Dec 28, 2023

Conversation

alexander-lakocy
Copy link

Description Of Changes

Adds calculation for Galvez-Davison Index for tropical convection.

Checklist

This is all visible in the changelist but wanted to document the things I had to change. I didn't see any of this in the Contributor's Guide, ideally we can add it under documentation to make it easier for contributors to document future feature additions.

@alexander-lakocy alexander-lakocy requested a review from a team as a code owner July 16, 2023 20:47
@alexander-lakocy alexander-lakocy requested review from dopplershift and removed request for a team July 16, 2023 20:47
@CLAassistant
Copy link

CLAassistant commented Jul 16, 2023

CLA assistant check
All committers have signed the CLA.

@dopplershift dopplershift added this to the September 2023 milestone Aug 30, 2023
@kgoebber
Copy link
Collaborator

Thanks for the contribution @alexander-lakocy! I was doing some quick validation and found the GDI calculated via NCEP and for 00 UTC on 14 September 2023 get:
Screenshot 2023-09-14 at 2 52 40 PM

Then using the calculation provided in this PR I get the following:

gdi_calculation_example

Generally the same look and shape to areas, but the magnitude appears to be a bit different. Do you know of someplace that has this calculated on a grid that would be available for a more detailed comparison (e.g., beyond the visual look at two different plots)?

@alexander-lakocy
Copy link
Author

Generally the same look and shape to areas, but the magnitude appears to be a bit different. Do you know of someplace that has this calculated on a grid that would be available for a more detailed comparison (e.g., beyond the visual look at two different plots)?

I agree the magnitude seems off. It would be useful to have a set of calculated values for creating a useful test. There were not any data grids in the source paper (although I believe they had some sample maps). Maybe we can check with the issue author to see if they have any example test problems?

@kgoebber
Copy link
Collaborator

kgoebber commented Sep 14, 2023

Okay, I think I found it...
gdi_calculation_example

It turns out that with the GFS it was using the p_sfc = 1 Pa!

The calculation calls for the surface pressure, so when calculating the value for a 1D vertical profile with monotonically increasing values, the current PR captures that with the line

p_sfc = pressure[0]

Which would be all good.

However, when calculating with the values from the GFS the pressure values are organized in the opposite way with the lowest pressure first and the highest pressure last. Additionally, the maximum value from that will not give the true surface pressure as that would be a separate variable. So we have a choice as to whether to mandate the input of the surface pressure explicitly (the most robust option), or have that be optional and assume a value by using the vertical pressure given and finding the maximum value (for a grid that would likely be 1000 hPa). I think I lean toward adding a fourth parameter to be true to the calculation and that would just add an extra input for someone doing the calculation on a 1D profile.

Additionally, the paper states that the gamma parameter should have units of 1/K, but it really should be 1/K^2, which them makes all the units work out and there would be no need to get the magnitude of all of the components and attach dimensional units at the end. I'll put a few comments in the code with some of the changes.

Copy link
Collaborator

@kgoebber kgoebber left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Lots of small requested changes that identify the issues brought up in a previous comment.

src/metpy/calc/thermo.py Outdated Show resolved Hide resolved
src/metpy/calc/thermo.py Outdated Show resolved Hide resolved
src/metpy/calc/thermo.py Outdated Show resolved Hide resolved
src/metpy/calc/thermo.py Outdated Show resolved Hide resolved
src/metpy/calc/thermo.py Outdated Show resolved Hide resolved
src/metpy/calc/thermo.py Outdated Show resolved Hide resolved
tests/calc/test_thermo.py Outdated Show resolved Hide resolved
tests/calc/test_thermo.py Show resolved Hide resolved
tests/calc/test_thermo.py Show resolved Hide resolved
src/metpy/calc/thermo.py Show resolved Hide resolved
@kgoebber
Copy link
Collaborator

@alexander-lakocy just a quick ping to see if you would like to update this PR with the requested changes or would you rather we take care of the requested changes?

@alexander-lakocy
Copy link
Author

alexander-lakocy commented Sep 21, 2023 via email

Copy link
Member

@dcamron dcamron left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With @kgoebber's additions, I'll step in to review the original code and the updates. A question to kick it off, I'll be back for more!

src/metpy/calc/thermo.py Show resolved Hide resolved
@kgoebber
Copy link
Collaborator

@dcamron reviewing the references within the GDI paper to dig into the value for C_p and L_v, it would appear most appropriate to use the constant values that are a part of MetPy. I don't imagine that there will be a very large difference in the calculated GDI, but the reference to Betts and Dugan (1974) and from there to Hilton (1972) would indicate the typical values that are traditionally used would be ~1004 J/kg/K and 2.5e6 J/kg for C_p and L_v, respectively.

@dcamron
Copy link
Member

dcamron commented Dec 15, 2023

but the reference to Betts and Dugan (1974) and from there to Hilton (1972) would indicate the typical values that are traditionally used would be ~1004 J/kg/K and 2.5e6 J/kg for C_p and L_v, respectively.

Perfect, thanks for digging that up!

@dopplershift
Copy link
Member

There's now a conflict here with test_thermo.py.

@kgoebber
Copy link
Collaborator

okay, I rebased and resolved the conflict.

@dopplershift
Copy link
Member

Tests are failing now, not sure why?

@kgoebber
Copy link
Collaborator

I think what we might be seeing there is the difference between the C_p values used. Currently this PR is updated to use the MetPy value and not the value from the paper. What may have happened is that I switched to using the MetPy value after I last pushed to the branch, but didn't actually follow through and update the tests.

@kgoebber
Copy link
Collaborator

I've updated the test values. I think the old values for the xarray test were the ones that used the wrong surface/starting pressure. I did a quick test with the different values of Cp_d and Lv. The difference with Cp_d are small, a few tenths of the GDI value. The difference with the choice of Lv is much more substantial. I believe the paper uses a value greater than 2.5e6 J/kg because of the index being used for a tropical atmosphere and with the temperature dependence on the latent heat of vaporization they chose not to use the value valid for 0C. The current implementation uses the value from the paper for Lv (which is the variable l_0 in the code/paper).

@dcamron
Copy link
Member

dcamron commented Dec 18, 2023

I am alright with that given the tropical context, and thanks for your work to investigate that. I'm going to take a small crack at some implementation changes today, but nothing to hold it up on if I don't quite get there. I'll ping Ryan for a final review then.

kgoebber and others added 4 commits December 26, 2023 12:01
Apply subjective code style cleanup to reduce repetition. Test 950
level check sooner. Collapse array dimension handling.
@dcamron dcamron force-pushed the feat/gdi branch 2 times, most recently from 0616bbb to 893a602 Compare December 26, 2023 19:06
@dopplershift dopplershift dismissed kgoebber’s stale review December 28, 2023 03:11

Updated and contributed to since this review.

@dopplershift dopplershift merged commit 2f8dc02 into Unidata:main Dec 28, 2023
34 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Please add GDI (Galvez Davison Index) for tropical convection
6 participants