Calculation of multi-level wet bulb temperature diagnostic

[bayliffe]

Following Tim's work on this ticket, we now need to implement two things I believe:

1. A method for generating a lookup table using the Goff-Gratch saturated vapour pressure equation. (A CLI perhaps?)
2. A Plugin that uses this lookup table to efficiently return the wet-bulb temperature from inputs of temperature, pressure and humidity. (Input cubes and output cube will have same dimensions)

I envisage the lookup table being a NetCDF file containing a cube with dimensions of temperature, pressure and humidity and the method making use of scipy.interpolate.interpnd.

====

Part of Snow melt Epic #191

As an IMPROVER developer, I need a multi-level wet bulb temperature diagnostic to be used in the calculation of the snow melting level.

Example formula from NOAA website. https://www.weather.gov/media/epz/wxcalc/wetBulbTdFromRh.pdf

NB the link above essentially requires the use of a piece of graph paper as far as I can see? As an alternative, please see UM documentation UMDP 080 Dynamical and Physical Diagnostic Calculations, page 4.

Acceptance criteria:

• Calculation of multi-level wet bulb temperature diagnostic from temperature, pressure, and relative humidity.
• Calculated on standard grid vertical levels.
• Need levels up to ~2000m for the UK.
• Compare values with UM model output given by stash 09222 (WET BULB TEMPERATURE (K)).
• Take a note of speed of calculation.
• Review of methodology

[gavinevans]

As part of being able to calculate the wet bulb temperature from temperature, pressure and relative humidity, the following intermediate diagnostics may be required:

• water vapour pressure
• saturated water vapour pressure
• mixing ratio
• saturated mixing ratio.

[PaulAbernethy]

I found this while doing a quick google search on the issue. Might provide a good starting point: https://github.com/remcmurry/Psychropy/blob/master/psychropy.py

[wxtim]

The story so far...

I have created a branch https://github.com/wxtim/improver/tree/IMPRO-174_psychrometric_calculations_wet_bulb and produced

method	function	unit tests
saturated_vapour_pressure	Yes	full set
humidity_ratio	Yes	basic functionality
wet_bulb_temperature	No	No
overall	No	No

I have a number of concerns about the code so far, and there is a good deal in the comments section describing things I'm not happy about. In brief:

• Is temp.convert_units('K') potentially going to slow this code down because it loads the whole cube into memory?
• Are we using an overly complex algorithm that we don't know enough about. What's wrong with the much simpler formula at https://www.weather.gov/media/epz/wxcalc/vaporPressure.pdf
• The algorithm is very unreadable atm. Some tidying up might be in order.
• At the moment we are pulling the data out of its Iris cube wrapper and operating on it. It is possible that we might prefer to operate on the cubes, although it's somewhat fiddly.
• Might our user wish to give pressure as a constant float rather than a cube?

[wxtim]

@PaulAbernethy has suggested that I look at papers:

Numerical data and functional relationships in science and technology. New series. Group V. Volume 4. Meteorology. Subvolume b. Physical and chemical properties of the air, P35 for the Goff-Gratch equation for saturated vapour pressure.

Gill, Atmosphere-Ocean Dynamics, Appendix 4 Equation A4.7: For the conversion of saturated vapour pressure over water to saturated vapour pressure in air. ScanToMyDocs_20170821.pdf

[wxtim]

State of Play:

I have created a pull request #261 for this ticket, although I don't think I am nearly there yet.

The basic structure of the calculations encapsulated in the pull request is shown in the attached picture. What this shows is that there is a web of dependency between different algorithms. Most of these are formulae based on regression using as many as 13 undistinguished constants.

Saturation Vapour Pressure

I have based wet bulb calculations upon the Goff-Gratch Saturated Vapour Pressure equation because this is how the WMO recommend carrying out this calculation. None of the methods is noticeably computationally intensive when calculating a few hundred elements for a comparison. The results of the simplistic method diverge significantly from the other two for high but meteorlogically possible temperatures, so I have included it only for completeness.

Wet Bulb Temperature

In pseudocode:

calculate the humidity mixing ratio using relative humidity
calculate the humidity mixing ratio using temperature = wetbulb temperature
do:
    calculate the humidity mixing ratio based on your new value of wet-bulb temperature
while: difference between your two estimates of humidity mixing ratio is greater than ...?

This iterative method will give answers in a group of extreme cases in the unit tests, but it displays a tendency to produce -inf and nans if used outside these ranges. I am still investigating what range of temperatures and pressures it should be used for.

Outstanding work

• [ ] These methods may require refactoring into a plugin method
• [ ] The saturated vapour pressure functions may require a wrapping function to allow the user to select the desired method through a kwarg.
• [x] I am considering how to make the unit tests more thorough, and how to catch issues.
• [x] Check how this handles 3D cubes.

[wxtim]

Summary of Acceptance Criteria

• [x] Calculation of multi-level wet bulb temperature diagnostic from temperature, pressure, and relative humidity. Results listed in jupyter notebook at Link on Jira ticket
• [x] Calculated on standard grid vertical levels - Need levels up to ~2000m for the UK. - Cannot easily test this with the data available in the sample data, but I don't see why it won't work, unit tests written to cover this
• [ ] Compare values with UM model output given by stash 09222 (WET BULB TEMPERATURE (K)). I'm not sure this data is available at the moment - I cannot find it looking through /critical model output.
• [x] Take a note of speed of calculation. - It's very slow -~ 530 seconds a cube of size (realization: 2; projection_y_coordinate: 970; projection_x_coordinate: 1042)
• [x] Review of methodology. - How formally do we want to write this up?

Summary of DOD

• [ ] Meets acceptance criteria in issue (if it is a bugfix PR, the bugfix PR should accurately spell out the problem being fixed with a user story). - Not yet
• [ ] Meets code style guide.
• [ ] All new plugin functionality unit tested if it is reasonable to do so (i.e. if the code is anything other than obvious).
• [x] All tests pass (should be run locally using /bin/improver-tests as well as checking travis passes).
• [x] Codacy GitHub check passes (if Codacy raises something unreasonable, ask for a tweak to be made to Codacy settings to allow it to pass)
• [ ] Known bugs fixed, or follow-on issues have been raised. Make sure the product owner is aware of any follow on tickets that have been raised.
• [ ] Reviewed by at least one person, but always at least 2 for serious non-bugfix changes.
• [x] Ensure licence information referred to in the code style guide is included within any new files.
• [x] Prepare input to sprint demo e.g. slides.

[wxtim]

[wxtim]

Review of Methodology.pdf

[wxtim]

It has been pointed out to me that the most efficient way of doing this is to use the code I have created to make an ancillary lookup table/cube with dimensions (temperature, relative_humidity, pressure).

This is very much analogous to how the UM code works.

the following module may be of use scipy.interpolate.interpnd

[MoseleyS]

Updated description to show what still needs to be done for this ticket.

[gavinevans]

Closing this issue, as this work will be carried on in #290 and #291, now that we've agreed an approach.