Potential Layer Selection Error in Storm Motion Calculation (`wrf_relhl.f90`)

[yfNJU]

Issue Description:​

In wrf_relhl.f90, the storm motion vector is estimated as 75% of the mean wind speed and 30° to the right of the mean wind direction (lines 103-113). This method is similar to Maddox et al. (1976), which computes the mean wind using surface, 850-, 700-, 500-, 300-, and 200-hPa levels. However, a potential layer selection issue exists in the mean wind calculation between 3000–10,000 m AGL:

1. The code identifies k3(layer at/above 3000 m AGL) and k10(layer at/above 10,000 m AGL) (lines 70–81)
2. The current loop DO k = k3, k10, -1(lines 87–92) likely includes data abovethe target heights, not strictly between3000–10,000 m AGL.

Proposed Fix:​​

Change the loop to DO k = k3+1, k10+1, -1 to ensure averaging only layers within 3000–10,000 m AGL

​Questions:​​

1. Why is the 3000–10,000 m AGL layer used? Is the current inclusion of heights outside this range intentional?
2. Is the proposed fix appropriate to align with the physical intent?

Thanks!

[kafitzgerald]

Thanks for the question.

I'm not sure offhand, but may be able to fill in some gaps at least.

Much of the Fortran code in here came from NCL, which is also on GitHub and it looks like this particular code was pulled from RIP4 into NCL according to the commit history.

The RIP codebase and relevant file can be found here: https://github.com/NCAR/RIP/blob/master/src/relhl.f
And the corresponding documentation here: https://www2.mmm.ucar.edu/wrf/users/docs/ripug.htm
Unfortunately, the commit history doesn't seem to have been retained in the RIP GitHub repo.

As for the inclusion of the layers above/below 3,000-10,000 in the averaging calculations, I suspect this was intentional in order to average over something close to the full range. Though as you pointed out, it does include values from outside of the range in the calculations.

[yfNJU]

Thank you for your reply! I believe modifying the k index to DO k = k3+1, k10+1, -1 would be more appropriate. May I submit a pull request for this? Additionally, I think using the Bunker's Internal Dynamics method to calculate storm motion is also a common approach. Perhaps we could also incorporate the relevant calculation method? Thanks!

[kafitzgerald]

Just to clarify, it looks like the current implementation averages values from both above and below 3000 and 10000 at either end of the do loop respectively. Presumably this is intended to approximate something near those values, so while it's not limited to including values only from within the specified range it doesn't seem incorrect. In order to average only over values from strictly within the 3000-10000 range, the indexing would need to look something like DO k = k3+1, k10-1, -1 and I'm not sure we'd want to do that especially given the long history of this algorithm being written as is.

Thanks for the recommendation regarding Bunker's method. It looks like this is something MetPy may have already implemented (i.e. bunkers_storm_motion). Is that something that might meet your needs?

Unfortunately, there's very limited support/bandwidth available for WRF-Python at this point. So while we would consider a pull request for new functionality, I'm not sure that implementing things here makes the most sense at this point and would suggest looking at other packages that integrate with the broader Xarray / Pangeo ecosystem (e.g. MetPy, GeoCAT packages, etc.).

That said, if you have requests for Earth system science / atmospheric science functionality that doesn't yet exist or meet your needs, we would be happy to take suggestions via GitHub or this feedback form.