improver.psychrometric_calculations.psychrometric_calculations module
Module to contain Psychrometric Calculations.
- class HumidityMixingRatio(model_id_attr=None)[source]
Bases:
BasePluginReturns the humidity mass mixing ratio from temperature, pressure and relative humidity
- _abc_impl = <_abc_data object>
- class PhaseChangeLevel(phase_change, grid_point_radius=2, horizontal_interpolation=True, model_id_attr=None)[source]
Bases:
BasePluginCalculate a continuous field of heights relative to sea level at which a phase change of precipitation is expected.
- __init__(phase_change, grid_point_radius=2, horizontal_interpolation=True, model_id_attr=None)[source]
Initialise class.
- Parameters:
phase_change (
str) –The desired phase change for which the altitude should be returned. Options are:
snow-sleet - the melting of snow to sleet. sleet-rain - the melting of sleet to rain. hail-rain - the melting of hail to rain.
grid_point_radius (
int) – The radius in grid points used to calculate the maximum height of the orography in a neighbourhood to determine points that should be excluded from interpolation for being too close to the orographic feature where high-resolution models can give highly localised results. Zero uses central point only (neighbourhood is disabled). One uses central point and one in each direction. Two goes two points etc. A grid_point_radius may be specified for data on any projection but the effective kernel shape in real space may be irregular. Users must be aware of this when choosing whether to use a non-zero grid_point_radius with a non-equal areas projectionhorizontal_interpolation (
bool) – If True apply horizontal interpolation to fill in holes in the returned phase-change-level that occur because the level falls below the orography. If False these areas will be masked.model_id_attr (str) – Name of the attribute used to identify the source model for blending.
- _abc_impl = <_abc_data object>
- _calculate_phase_change_level(wet_bulb_temp, wb_integral, orography, max_nbhood_orog, land_sea_data, heights, height_points, highest_height)[source]
Calculate phase change level and fill in missing points
Method
The first call is to
find_falling_levelwhich finds the height at which the wet-bulb integral reaches the required threshold, this threshold defining the phase change level. The threshold value depends upon the phase change that has been requested. If we are unable to find the height of the phase level for a grid point from the wet-bulb integral the value at that point is set to np.nan and will be filled in later.The next call is to
fill_in_high_phase_change_falling_levels. This function is used to fill in any data in the phase change level where this level is very high, effectively above the highest level used in the wet bulb temperature integration. In these areas the atmosphere is so warm that the phase change has occured above the highest level that is being considered. In these areas the phase change level is set to the highest height level + the height of the orography. The exact height of the phase change level is unimportant in these areas as we can be certain the precipitation at the surface will be rain.Next we fill in any sea points where we have not found a phase change level by the time we get to sea level, i.e. where the wet bulb temperature integral never reaches the desired threshold. To fill these points we call
fill_in_sea_pointswhich finds a linear fit to the wet bulb temperature close to sea level and uses this to find where an extrapolated wet bulb temperature integral would cross the threshold. This results in phase change levels below sea level for points where we have applied the extrapolation. It is important to get sensible values for these heights, particularly in marginal snow fall cases. As we calculate probabilities using multiple realizations, if some realizations give a snow falling level just above sea level and some just below, getting these heights close to correct ensures we calculate realistic probabilities of snow falling at the surface. Consider instead if we simply set the falling levels at these points to be equal to sea level; in so doing we would be introducing a warm bias to the output.Finally, if there are any areas for which the phase change level remains unset, we know that these must be at or below the height of the orography. As such we fill these in using horizontal interpolation across the grid. This is achieved using two calls to
interpolate_missing_data. The first attempts to use linear interpolation to fill as many of these holes as possible. The second uses nearest neighbour filling for any areas where linear interpolation was not possible, e.g. at the edge of the domain. The process is as follows:Fill in the phase change level for points with no value yet set using horizontal interpolation from surrounding set points. Only interpolate from surrounding set points at which the phase change level is below the maximum orography height in the region around the unset point. This helps us avoid spreading very high phase change levels across areas where we had missing data.
Fill any gaps that still remain where the linear interpolation has not been able to find a value because there is not enough data (e.g at the corners of the domain). Use nearest neighbour filling.
Check whether despite our efforts we have still filled in some of the missing points with phase change levels above the orography. In these cases set the missing points to the height of orography.
If there are any points that have somehow remain unset, these are filled with the missing_data value.
- Parameters:
wet_bulb_temp (
ndarray) – Wet bulb temperature datawb_integral (
ndarray) – Wet bulb temperature integralorography (
ndarray) – Orography heightsmax_nbhood_orog (
ndarray) – Maximum orography height in neighbourhood (used to determine points that can be used for interpolation)land_sea_data (
ndarray) – Mask of binary land / sea dataheights (
ndarray) – All heights of wet bulb temperature inputheight_points (
ndarray) – Heights on wet bulb temperature integral slicehighest_height (
float) – Height of the highest level to which the wet bulb temperature has been integrated
- Return type:
- Returns:
Level at which phase changes
- _horizontally_interpolate_phase(phase_change_data, orography, max_nbhood_orog)[source]
Fill in missing points via horizontal interpolation.
- Parameters:
- Return type:
- Returns:
Level at which phase changes, with missing data filled in
- create_phase_change_level_cube(wbt, phase_change_level)[source]
Populate output cube with phase change data
- fill_in_high_phase_change_falling_levels(phase_change_level_data, orog_data, highest_wb_int_data, highest_height)[source]
Fill in any data in the phase change level where the whole wet bulb temperature integral is above the the threshold. Set these points to the highest height level + orography.
- fill_in_sea_points(phase_change_level_data, land_sea_data, max_wb_integral, wet_bulb_temperature, heights, orography)[source]
Fill in any sea points where we have not found a phase change level by the time we get to sea level, i.e. where the whole wet bulb temperature integral is below the threshold.
This function finds a linear fit to the wet bulb temperature close to sea level and uses this to find where an extrapolated wet bulb temperature integral would cross the threshold. This results in phase change levels below sea level for points where we have applied the extrapolation.
The linear fit assumes that all sea points have an orography altitude of zero, however this is not always the case. The orography altitude is added to phase change levels calculated from the linear fit to account for this.
Assumes that height is the first axis in the wet_bulb_integral array.
- Parameters:
phase_change_level_data (
ndarray) – The phase change level array, filled with values for points whose wet bulb temperature integral crossed the theshold.land_sea_data (
ndarray) – The binary land-sea maskmax_wb_integral (
ndarray) – The wet bulb temperature integral at the final height level used in the integration. This has the maximum values for the wet bulb temperature integral at any level.wet_bulb_temperature (
ndarray) – The wet bulb temperature profile at each grid point, with height as the leading dimension.heights (
ndarray) – The vertical height levels above orography, matching the leading dimension of the wet_bulb_temperature.orography (
ndarray) – Orography heights
- Return type:
- find_extrapolated_falling_level(max_wb_integral, gradient, intercept, phase_change_level_data, sea_points)[source]
Find the phase change level below sea level using the linear extrapolation of the wet bulb temperature integral and update the phase change level array with these values.
The phase change level is calculated from finding the point where the integral of wet bulb temperature crosses the falling level threshold.
In cases where the wet bulb temperature integral has not reached the threshold by the time we reach sea level, we can find a fit to the wet bulb temperature profile near the surface, and use this to estimate where the phase change level would be below sea level.
The difference between the wet bulb temperature integral at the threshold and the wet bulb integral at the surface is equal to the integral of the wet bulb temperature between sea level and the negative height corresponding to the phase change level. As we are using a simple linear fit, we can integrate this to find an expression for the extrapolated phase change level.
The form of this expression depends on whether the linear fit of wet bulb temperature crosses the height axis above or below zero altitude.
If we have our linear fit of the form:
\[{{wet\:bulb\:temperature} = m \times height + c}\]and let \(I\) be the wet bulb temperature integral we have found above sea level.
If it crosses above zero, then the limits on the integral are the phase change level and zero and we find the following expression for the phase change level:
\[{{phase\:change\:level} = \frac{c \pm \sqrt{ c^2-2 m (threshold-I)}}{-m}}\]If the linear fit crosses below zero the limits on our integral are the phase change level and the point where the linear fit crosses the height axis, as only positive wet bulb temperatures count towards the integral. In this case our expression for the phase change level is:
\[{{phase\:change\:level} = \frac{c \pm \sqrt{ 2 m (I-threshold)}}{-m}}\]- Parameters:
max_wb_integral (
ndarray) – The wet bulb temperature integral at sea level.gradient (
ndarray) – The gradient of the line of best fit we are using in the extrapolation.intercept (
ndarray) – The intercept of the line of best fit we are using in the extrapolation.phase_change_level_data (
ndarray) – The phase change level array with values filled in with phase change levels calculated through extrapolation.sea_points (
ndarray) – A boolean array with True where the points are sea points.
- Return type:
- find_falling_level(wb_int_data, orog_data, height_points)[source]
Find the phase change level by finding the level of the wet-bulb integral data at the required threshold. Wet-bulb integral data is only available above ground level and there may be an insufficient number of levels in the input data, in which case the required threshold may lie outside the Wet-bulb integral data and the value at that point will be set to np.nan.
- find_max_in_nbhood_orography(orography_cube)[source]
Find the maximum value of the orography in the neighbourhood around each grid point. If self.grid_point_radius is zero, the orography is used without neighbourhooding.
- static linear_wet_bulb_fit(wet_bulb_temperature, heights, sea_points, start_point=0, end_point=5)[source]
Calculates a linear fit to the wet bulb temperature profile close to the surface to use when we extrapolate the wet bulb temperature below sea level for sea points.
We only use a set number of points close to the surface for this fit, specified by a start_point and end_point.
- Parameters:
wet_bulb_temperature (
ndarray) – The wet bulb temperature profile at each grid point, with height as the leading dimension.heights (
ndarray) – The vertical height levels above orography, matching the leading dimension of the wet_bulb_temperature.sea_points (
ndarray) – A boolean array with True where the points are sea points.start_point (
int) – The index of the the starting height we want to use in our linear fit.end_point (
int) – The index of the the end height we want to use in our linear fit.
- Return type:
- Returns:
An array, the same shape as a 2D slice of the wet_bulb_temperature input, containing the gradients of the fitted straight line at each point where it could be found, filled with zeros elsewhere.
An array, the same shape as a 2D slice of the wet_bulb_temperature input, containing the intercepts of the fitted straight line at each point where it could be found, filled with zeros elsewhere.
- process(cubes)[source]
Use the wet bulb temperature integral to find the altitude at which a phase change occurs (e.g. snow to sleet). This is achieved by finding the height above sea level at which the integral matches an empirical threshold that is expected to correspond with the phase change. This empirical threshold is the falling_level_threshold. Fill in missing data appropriately.
- Parameters:
containing (cubes) –
- wet_bulb_temperature:
Cube of wet bulb temperatures on height levels.
- wet_bulb_integral:
Cube of wet bulb temperature integral (Kelvin-metres).
- orog:
Cube of orography (m).
- land_sea_mask:
Cube containing a binary land-sea mask, with land points set to one and sea points set to zero.
- Return type:
- Returns:
Cube of phase change level above sea level (asl).
- Raises:
ValueError – Raise exception if the model_id_attr attribute does not match on the input cubes.
- _calculate_latent_heat(temperature)[source]
Calculate a temperature adjusted latent heat of condensation for water vapour using the relationship employed by the UM.
- _latent_heat_release(q1, q2, temperature)[source]
Returns the latent heat released (K) when condensing water vapour from specific humidity value q1 to q2, both in kg kg-1 when the temperature is approximately t. Returns negative values when initial condition is subsaturated as this method assumes there is always liquid water present which can be evaporated.
- Parameters:
- Return type:
- Returns:
Temperature adjustment to apply to account for latent heat release (K).
- _svp_from_lookup(temperature)[source]
Gets value for saturation vapour pressure in a pure water vapour system from a pre-calculated lookup table. Interpolates linearly between points in the table to the temperatures required.
- _svp_table()[source]
Calculate a saturated vapour pressure (SVP) lookup table. The lru_cache decorator caches this table on first call to this function, so that the table does not need to be re-calculated if used multiple times.
A value of SVP for any temperature between T_MIN and T_MAX (inclusive) can be obtained by interpolating through the table, as is done in the _svp_from_lookup function.
- Return type:
- Returns:
Array of saturated vapour pressures (Pa).
- adjust_for_latent_heat(temperature_in, humidity_in, pressure)[source]
Increases temperature and reduces humidity via latent heat release from condensation until values represent 100% relative humidity.
Subsaturated values will be returned unaltered.
This method uses the scipy newton solver with a limit of 6 iterations. The deepest convection needs more iterations to converge. This is only important if we reach the position that all points in an array fail to converge at the same pressure level, because the solver raises an exception (although docs say it shouldn’t). I haven’t seen this be true except when the array contains only one point, but increasing the maximum number of iterations for small arrays is a very small price to pay for stability.
The latent heat energy \(Q\) released through condensation of water is equal to the mass of water condensed (\(m\)) multiplied by the constant latent heat of condensation of water (\(L\)) defined in the improver.constants module.
\[Q = L m\]The temperature change of the air parcel is the latent heat energy divided by the specific heat of dry air (\(C_p\)) defined in the improver.constants module.
\[\Delta T = \frac{Q}{C_p}\]The mass of water condensed is dependent on the difference between the initial parcel temperature and the final parcel temperature, because the saturated mass mixing ratio \(q_{sat}\) of water in air is dependent on the temperature and pressure of the air. The derivation of \(q_{sat}\) is in
saturated_humidity().Therefore, when ascending a saturated parcel of air from pressure level \(P_1\) to \(P_2\), the resulting temperature is
\[T_2 = T_1 + \frac{L}{C_p}(q_{sat(T_1,P_1)} - q_{sat(T_2,P_2)})\]This requires an iterative solver.
The solver takes in the initial value of the humidity mixing ratio (\(q_1\)), the dry-adiabatically adjusted pressure (\(P_2\)) and temperature (\(T_2\)) (see
dry_adiabatic_temperature()), and a guess of the final saturated humidity mixing ratio (\(q_{sat(T_{2}^{'})}\)). From these, it calculates the mass of water unaccounted for:\[\Delta q = q_{sat(T_2,P_2)} - q_{sat(T_{2}^{'},P_2)}\]Where
\[T_{2}^{'} = T_2 + \frac{L}{C_p}(q_{sat(T_2,P_2)} - q_{sat(T_{2}^{'},P_2)})\]The solver iterates until \(\Delta q < 1 \times 10^{-6}\). The final temperature (\(T_{2}^{'}\)) is calculated from the
saturated_humidity()method and returned along with the final humidity mixing ratio.- Parameters:
- Return type:
- Returns:
tuple of temperature (K) and humidity (kg kg-1) after saturated latent heat release
- calculate_svp_in_air(temperature, pressure)[source]
Calculates the saturation vapour pressure in air. Looks up the saturation vapour pressure in a pure water vapour system, and pressure-corrects the result to obtain the saturation vapour pressure in air.
- Parameters:
- Return type:
- Returns:
Saturation vapour pressure in air (Pa).
References
Atmosphere-Ocean Dynamics, Adrian E. Gill, International Geophysics Series, Vol. 30; Equation A4.7.
- dry_adiabatic_pressure(initial_temperature, initial_pressure, final_temperature)[source]
Calculate pressure at final_temperature after adiabatic adjustment of dry air from the initial temperature and pressure.
The dry adiabatic equation for finding the target pressure is
\[ \begin{align}\begin{aligned}P_2 = P_1 \left(\frac{T_2}{T_1}\right)^\frac{1}{\kappa}\\\kappa = \frac{R}{C_p}\end{aligned}\end{align} \]Where \(R\) is the gas constant and \(C_p\) is the specific heat of dry air. Both are defined in the improver.constants module
- dry_adiabatic_temperature(initial_temperature, initial_pressure, final_pressure)[source]
Calculate temperature at final_pressure after adiabatic adjustment of dry air from the initial temperature and pressure.
The dry adiabatic equation for finding the target temperature is
\[ \begin{align}\begin{aligned}T_2 = T_1 \left(\frac{P_2}{P_1}\right)^\kappa\\\kappa = \frac{R}{C_p}\end{aligned}\end{align} \]Where \(R\) is the gas constant and \(C_p\) is the specific heat of dry air. Both are defined in the improver.constants module
- saturated_humidity(temperature, pressure)[source]
Calculate specific humidity mixing ratio of saturated air of given temperature and pressure
- Parameters:
- Return type:
- Returns:
Array of specific humidity values (kg kg-1) representing saturated air
Method from referenced documentation. Note that EARTH_REPSILON is simply given as an unnamed constant in the reference (0.62198).
References
ASHRAE Fundamentals handbook (2005) Equation 22, 24, p6.8