UFO
Formulas.cc
Go to the documentation of this file.
1 /*
2  * (C) Crown copyright 2020, Met Office
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  */
7 
8 #include <algorithm>
9 #include <cmath>
10 #include <string>
11 #include <vector>
12 #include "eckit/exception/Exceptions.h"
13 #include "oops/util/Logger.h"
14 #include "ufo/utils/Constants.h"
17 
18 namespace ufo {
19 
20 namespace formulas {
21 
22 MethodFormulation resolveMethods(const std::string& input) {
23  // Met Centers
24  if (input == "UKMO") return UKMO;
25  if (input == "NCAR") return NCAR;
26  if (input == "NOAA") return NOAA;
27  return DEFAULT;
28 }
29 
30 MethodFormulation resolveFormulations(const std::string& input, const std::string& method) {
31  if (input.empty()) {
32  // if we do not have any formulation set, then we assign
33  // the method as formulation
34  return resolveMethods(method);
35  } else {
36  // Formulation authors
37  if (input == "Sonntag") return Sonntag;
38  if (input == "Walko") return Walko;
39  if (input == "Murphy") return Murphy;
40  return DEFAULT;
41  }
42 }
43 
44 /* -------------------------------------------------------------------------------------*/
45 float SatVaporPres_fromTemp(float temp_K, MethodFormulation formulation) {
46  const float missingValueFloat = util::missingValue(1.0f);
47  const float t0c = static_cast<float>(ufo::Constants::t0c);
48  float e_sub_s = missingValueFloat; // Saturation vapour pressure (Pa)
49 
50  switch (formulation) {
53  /* I. Source: Eqn 7, Sonntag, D., Advancements in the field of hygrometry,
54  * Meteorol. Zeitschrift, N. F., 3, 51-66, 1994.
55  * Most radiosonde manufacturers use Wexler, or Hyland and Wexler
56  * or Sonntag formulations, which are all very similar (Holger Vomel,
57  * pers. comm., 2011)
58  */
59  if (temp_K != missingValueFloat) {
60  e_sub_s = std::exp(-6096.9385f / temp_K + 21.2409642f - 2.711193E-2f * temp_K +
61  1.673952E-5f * temp_K * temp_K + 2.433502f * std::log(temp_K));
62  } else {
63  e_sub_s = 0.0f;
64  }
65  break;
66  }
68  /* Returns a saturation mixing ratio given a temperature and pressure
69  using saturation vapour pressures caluclated using the Goff-Gratch
70  formulae, adopted by the WMO as taken from Landolt-Bornstein, 1987
71  Numerical Data and Functional relationships in Science and
72  Technology. Group V/Vol 4B Meteorology. Physical and Chemical
73  properties of Air, P35.
74  */
75  if (temp_K != missingValueFloat) {
76  float adj_Temp;
77  float lookup_a;
78  int lookup_i;
79 
80  const float Low_temp_thd = 183.15; // Lowest temperature for which look-up table is valid
81  const float High_temp_thd = 338.15; // Highest temperature for which look-up table is valid
82  const float Delta_Temp = 0.1; // Temperature increment of look-up table
83 
84  // Use the lookup table to find saturated vapour pressure.
85  adj_Temp = std::max(Low_temp_thd, temp_K);
86  adj_Temp = std::min(High_temp_thd, adj_Temp);
87 
88  lookup_a = (adj_Temp - Low_temp_thd + Delta_Temp) / Delta_Temp;
89  lookup_i = static_cast<int>(lookup_a);
90  lookup_a = lookup_a - lookup_i;
91  e_sub_s = (1.0 - lookup_a) *
93  lookup_a *
95  } else {
96  e_sub_s = 0.0f;
97  }
98  break;
99  }
101  // Polynomial fit of Goff-Gratch (1946) formulation. (Walko, 1991)
102  float x = std::max(-80.0f, temp_K-t0c);
103  const std::vector<float> c{610.5851f, 44.40316f, 1.430341f, 0.2641412e-1f,
104  0.2995057e-3f, 0.2031998e-5f, 0.6936113e-8f, 0.2564861e-11f, -0.3704404e-13f};
105  e_sub_s = c[0]+x*(c[1]+x*(c[2]+x*(c[3]+x*(c[4]+x*(c[5]+x*(c[6]+x*(c[7]+x*c[8])))))));
106  break;
107  }
109  // ALTERNATIVE (costs more CPU, more accurate than Walko, 1991)
110  // Source: Murphy and Koop, Review of the vapour pressure of ice and
111  // supercooled water for atmospheric applications, Q. J. R.
112  // Meteorol. Soc (2005), 131, pp. 1539-1565.
113  e_sub_s = std::exp(54.842763f - 6763.22f / temp_K - 4.210f * std::log(temp_K)
114  + 0.000367f * temp_K + std::tanh(0.0415f * (temp_K - 218.8f))
115  * (53.878f - 1331.22f / temp_K - 9.44523f * std::log(temp_K)
116  + 0.014025f * temp_K));
117  break;
118  }
122  default: {
123  // Classical formula from Rogers and Yau (1989; Eq2.17)
124  e_sub_s = 1000. * 0.6112 * std::exp(17.67f * (temp_K - t0c) / (temp_K - 29.65f));
125  break;
126  }
127  }
128  return e_sub_s;
129 }
130 
131 /* -------------------------------------------------------------------------------------*/
132 float SatVaporPres_correction(float e_sub_s, float temp_K, float pressure,
133  MethodFormulation formulation) {
134  const float t0c = static_cast<float>(ufo::Constants::t0c);
135 
136  switch (formulation) {
141  /* e_sub_s is the saturation vapour pressure of pure water vapour.FsubW (~ 1.005
142  at 1000 hPa) is the enhancement factor needed for moist air.
143  If P is set to -1: then eg eqns 20, 22 of Sonntag is used
144  If P > 0 then eqn A4.6 of Adrian Gill's book is used to guarantee consistency with the
145  saturated specific humidity.
146  */
147  float FsubW; // Enhancement factor
148  FsubW = 1.0f + 1.0E-8f * pressure * (4.5f + 6.0E-4f * (temp_K - t0c) *(temp_K - t0c));
149  e_sub_s = e_sub_s * FsubW;
150 
151  break;
152  }
153  default: {
154  std::string errString = "Aborting, no method matches enum formulas::MethodFormulation";
155  oops::Log::error() << errString;
156  throw eckit::BadValue(errString);
157  }
158  }
159  return e_sub_s;
160 }
161 /* -------------------------------------------------------------------------------------*/
162 
163 float Qsat_From_Psat(float Psat, float P, MethodFormulation formulation) {
164  float QSat = util::missingValue(1.0f); // Saturated specific humidity or
165  // saturated vapour pressure (if P<0)
166 
167  switch (formulation) {
171  default: {
172  // Calculation using the Sonntag (1994) formula. (With fix at low
173  // pressure)
174  // Note that at very low pressures we apply a fix, to prevent a
175  // singularity (Qsat tends to 1.0 kg/kg).
176  QSat = (Constants::epsilon * Psat) /
177  (std::max(P, Psat) - (1.0f - Constants::epsilon) * Psat);
178  break;
179  }
180  }
181  return QSat;
182 }
183 
184 /* -------------------------------------------------------------------------------------*/
185 
186 // VirtualTemperature()
187 float VirtualTemp_From_Psat_P_T(float Psat, float P, float T, MethodFormulation formulation) {
188  float Tv = util::missingValue(1.0f); // virtual temperature
189 
190  switch (formulation) {
194  default: {
195  Tv = T * ((P + Psat / Constants::epsilon) / (P + Psat));
196  break;
197  }
198  }
199  return Tv;
200 }
201 
202 /* -------------------------------------------------------------------------------------*/
203 
204 float VirtualTemp_From_Rh_Psat_P_T(float Rh, float Psat, float P, float T,
205  MethodFormulation formulation) {
206  float Tv = util::missingValue(1.0f); // virtual temperature
207 
208  switch (formulation) {
212  default: {
213  float PsatRh = Psat * Rh * 0.01f;
214  Tv = VirtualTemp_From_Psat_P_T(PsatRh , P, T, formulation);
215  break;
216  }
217  }
218 
219  return Tv;
220 }
221 
222 /* -------------------------------------------------------------------------------------*/
223 
224 float Height_To_Pressure_ICAO_atmos(float height, MethodFormulation formulation) {
225  const float missingValueFloat = util::missingValue(1.0f);
226  float Pressure = missingValueFloat;
227 
228  switch (formulation) {
232  default: {
233  float RepT_Bot, RepT_Top, ZP1, ZP2;
234 
235  RepT_Bot = 1.0 / Constants::icao_temp_surface;
236  RepT_Top = 1.0 / Constants::icao_temp_isothermal_layer;
239 
240  if (height <= missingValueFloat) {
241  Pressure = missingValueFloat;
242  } else if (height < -5000.0) {
243  // TODO(david simonin): The original code has this test.
244  // Not sure why! Are we expecting very negative height value??
245  Pressure = missingValueFloat;
246  } else if (height < Constants::icao_height_l) {
247  // Heights up to 11,000 geopotential heigh in meter [gpm]
248  Pressure = Constants::icao_lapse_rate_l * height * RepT_Bot;
249  Pressure = std::pow((1.0 - Pressure), ZP1);
250  Pressure = 100.0 * Pressure * Constants::icao_pressure_surface;
251  } else if (height < Constants::icao_height_u) {
252  // Heights between 11,000 and 20,000 geopotential heigh in meter [gpm]
253  Pressure = Constants::g_over_rd * (height - Constants::icao_height_l) * RepT_Top;
254  Pressure = std::log(Constants::icao_pressure_l) - Pressure;
255  Pressure = 100.0 * std::exp(Pressure);
256  } else {
257  // Heights above 20,000 geopotential heigh in meter [gpm]
258  Pressure = Constants::icao_lapse_rate_u * RepT_Top *
259  (height - Constants::icao_height_u);
260  Pressure = 100.0 * Constants::icao_pressure_u *
261  std::pow((1.0 - Pressure), ZP2);
262  }
263  break;
264  }
265  }
266  return Pressure;
267 }
268 
269 /* -------------------------------------------------------------------------------------*/
270 
271 float Pressure_To_Height(float pressure, MethodFormulation method) {
272  const float missingValueFloat = util::missingValue(1.0f);
273  const float pressure_hPa = pressure * 0.01;
274  float height = missingValueFloat;
275 
276  switch (method) {
278  // The NCAR-RAL method: a fast approximation for pressures > 120 hPa.
279  // Above 120hPa (~15km) use the ICAO atmosphere.
280  if (pressure == missingValueFloat || pressure <= 0.0f) {
281  height = missingValueFloat;
282  } else if (pressure_hPa <= 120.0f &&
283  pressure_hPa > Constants::icao_pressure_u) {
284  pressure = std::log(Constants::icao_pressure_l) - std::log(pressure_hPa);
287  } else if (pressure_hPa <= Constants::icao_pressure_u) {
288  pressure = 1.0 - std::pow(pressure_hPa / Constants::icao_pressure_u,
292  } else {
293  height = 44307.692 * (1.0 - std::pow(pressure / 101325.0, 0.190));
294  }
295  break;
296 
299  default: {
300  if (pressure == missingValueFloat || pressure <= 0.0f) {
301  height = missingValueFloat;
302  } else if (pressure_hPa > Constants::icao_pressure_l) {
303  pressure = 1.0 - std::pow(pressure_hPa / Constants::icao_pressure_surface,
306  } else if (pressure_hPa <= Constants::icao_pressure_l &&
307  pressure_hPa > Constants::icao_pressure_u) {
308  pressure = std::log(Constants::icao_pressure_l) - std::log(pressure_hPa);
311  } else {
312  pressure = 1.0 - std::pow(pressure_hPa / Constants::icao_pressure_u,
316  }
317  break;
318  }
319  }
320  return height;
321 }
322 
323 float GetWindDirection(float u, float v) {
324  const float missing = util::missingValue(1.0f);
325  float windDirection = missing; // wind direction
326 
327  if (u != missing && v != missing) {
328  if (u != 0.0f || v != 0.0f) {
329  windDirection = fmod((270.0f - atan2(v, u) * Constants::rad2deg), 360.0f);
330  } else {
331  windDirection = 0.0f;
332  }
333  }
334  return windDirection;
335 }
336 
337 float GetWindSpeed(float u, float v) {
338  const float missing = util::missingValue(1.0f);
339  float windSpeed = missing; // wind speed
340 
341  if (u != missing && v != missing) {
342  windSpeed = hypot(u, v);
343  }
344  return windSpeed;
345 }
346 
347 float GetWind_U(float windSpeed, float windFromDirection) {
348  const float missing = util::missingValue(1.0f);
349  float u = missing; // wind speed
350 
351  if (windFromDirection != missing
352  && windSpeed != missing && windSpeed >= 0) {
353  u = -windSpeed * sin(windFromDirection * Constants::deg2rad);
354  }
355  return u;
356 }
357 
358 
359 float GetWind_V(float windSpeed, float windFromDirection) {
360  const float missing = util::missingValue(1.0f);
361  float v = missing; // wind speed
362 
363  if (windFromDirection != missing
364  && windSpeed != missing && windSpeed >= 0) {
365  v = -windSpeed * cos(windFromDirection * Constants::deg2rad);
366  }
367  return v;
368 }
369 
370 /* -------------------------------------------------------------------------------------*/
371 
372 int RenumberScanPosition(int scanpos) {
373  // Renumber from 2,5,8,... to 1,2,3,...
374  int newpos = (scanpos + 1)/3;
375  return newpos;
376 }
377 
378 /* -------------------------------------------------------------------------------------*/
379 
381 (const std::vector<size_t> & locs,
382  const std::vector<bool> & apply,
383  const std::vector<float> & lat_in,
384  const std::vector<float> & lon_in,
385  const std::vector<util::DateTime> & time_in,
386  const std::vector<float> & height,
387  const std::vector<float> & windspd,
388  const std::vector<float> & winddir,
389  std::vector<float> & lat_out,
390  std::vector<float> & lon_out,
391  std::vector<util::DateTime> & time_out,
392  MethodFormulation formulation) {
393  const float missingValueFloat = util::missingValue(1.0f);
394 
395  switch (formulation) {
399  default: {
400  // Location of the first entry in the profile.
401  const size_t loc0 = locs.front();
402 
403  // Values of latitude, longitude and datetime at the first entry of the profile.
404  const double lat0 = lat_in[loc0];
405  const double lon0 = lon_in[loc0];
406  const util::DateTime time0 = time_in[loc0];
407 
408  // The drift computation is not performed for very high latitude sites.
409  if (std::abs(lat0) >= 89.0) return;
410 
411  // Fill vector of valid locations.
412  std::vector<size_t> locs_valid;
413  for (size_t jloc : locs) {
414  // If not selected by the where clause.
415  if (!apply[jloc]) continue;
416  // The location is classed as valid if the wind speed and height are not missing.
417  if (windspd[jloc] != missingValueFloat && height[jloc] != missingValueFloat)
418  locs_valid.push_back(jloc);
419  }
420 
421  // If there are zero or one valid locations, exit the routine.
422  if (locs_valid.size() < 2) return;
423 
424  // Average ascent speed (m/s).
425  const double ascent_speed = 5.16;
426 
427  // Cumulative values of change in time.
428  // This value is converted to a util::Duration object rather than performing the
429  // same conversion to each individual change in time.
430  // This avoids a loss in precision given util::Duration is accurate to the nearest second.
431  double dt_cumul = 0.0;
432 
433  for (size_t k = 0; k < locs_valid.size() - 1; ++k) {
434  // Locations of the current and next valid observations in the profile.
435  const size_t loc_current = locs_valid[k];
436  const size_t loc_next = locs_valid[k + 1];
437 
438  // Compute changes in latitude, longitude and time between adjacent valid levels.
439  // Change in height.
440  const double dh = height[loc_next] - height[loc_current];
441  // Change in time.
442  const double dt = dh / ascent_speed;
443  // Average eastward and northward wind between the two levels.
444  // 180 degrees is subtracted from the wind direction in order to account for the different
445  // conventions used in the observations and in this calculation.
446  const double avgu = 0.5 *
447  (windspd[loc_current] * std::sin((winddir[loc_current] - 180.0) * Constants::deg2rad) +
448  windspd[loc_next] * std::sin((winddir[loc_next] - 180.0) * Constants::deg2rad));
449  const double avgv = 0.5 *
450  (windspd[loc_current] * std::cos((winddir[loc_current] - 180.0) * Constants::deg2rad) +
451  windspd[loc_next] * std::cos((winddir[loc_next] - 180.0) * Constants::deg2rad));
452  // Total height of the observation above the centre of the Earth.
453  const double totalheight = ufo::Constants::mean_earth_rad * 1000.0 + height[loc_current];
454  // Change in latitude.
455  const double dlat = ufo::Constants::rad2deg * avgv * dt / totalheight;
456  // Change in longitude.
457  const double dlon = ufo::Constants::rad2deg * avgu * dt /
458  (totalheight * std::cos(lat_out[loc_current] * ufo::Constants::Constants::deg2rad));
459 
460  // Fill output values.
461  lat_out[loc_next] = lat_out[loc_current] + dlat;
462  lon_out[loc_next] = lon_out[loc_current] + dlon;
463  // Convert the cumulative change in time to a util::Duration.
464  dt_cumul += dt;
465  time_out[loc_next] = time0 + util::Duration(static_cast<int64_t>(dt_cumul));
466  }
467 
468  // Copy latitude, longitude and time at each valid location to all invalid
469  // locations that lie between the current valid location and the next one above it.
470  double lat = lat0;
471  double lon = lon0;
472  util::DateTime time = time0;
473  for (size_t jloc : locs) {
474  if (std::find(locs_valid.begin(), locs_valid.end(), jloc) != locs_valid.end()) {
475  lat = lat_out[jloc];
476  lon = lon_out[jloc];
477  time = time_out[jloc];
478  } else {
479  lat_out[jloc] = lat;
480  lon_out[jloc] = lon;
481  time_out[jloc] = time;
482  }
483  }
484  break;
485  }
486  }
487 }
488 } // namespace formulas
489 } // namespace ufo
constexpr int missing
Definition: QCflags.h:20
float GetWind_U(float windSpeed, float windFromDirection)
Get eastward (u) wind component from wind speed and direction.
Definition: Formulas.cc:347
MethodFormulation resolveFormulations(const std::string &input, const std::string &method)
Definition: Formulas.cc:30
void horizontalDrift(const std::vector< size_t > &locs, const std::vector< bool > &apply, const std::vector< float > &lat_in, const std::vector< float > &lon_in, const std::vector< util::DateTime > &time_in, const std::vector< float > &height, const std::vector< float > &windspd, const std::vector< float > &winddir, std::vector< float > &lat_out, std::vector< float > &lon_out, std::vector< util::DateTime > &time_out, MethodFormulation formulation)
Compute horizontal drift latitude, longitude and time for an atmospheric profile.
Definition: Formulas.cc:381
float Pressure_To_Height(float pressure, MethodFormulation method)
Converts pressure to height.
Definition: Formulas.cc:271
float Qsat_From_Psat(float Psat, float P, MethodFormulation formulation)
Calculates Saturated specific humidity or saturated vapour pressure using saturation vapour pressure.
Definition: Formulas.cc:163
float VirtualTemp_From_Rh_Psat_P_T(float Rh, float Psat, float P, float T, MethodFormulation formulation)
Derive Virtual Tempreture using Relative humidity, sat.
Definition: Formulas.cc:204
float GetWindDirection(float u, float v)
Converts u and v wind component into wind direction.
Definition: Formulas.cc:323
float SatVaporPres_correction(float e_sub_s, float temp_K, float pressure, MethodFormulation formulation)
Calculates saturated vapour pressure from temperature.
Definition: Formulas.cc:132
float VirtualTemp_From_Psat_P_T(float Psat, float P, float T, MethodFormulation formulation)
Derive Virtual Temperature from saturation vapour pressure, pressure and temperature.
Definition: Formulas.cc:187
MethodFormulation resolveMethods(const std::string &input)
Definition: Formulas.cc:22
float SatVaporPres_fromTemp(float temp_K, MethodFormulation formulation)
Calculates saturated vapour pressure from temperature.
Definition: Formulas.cc:45
MethodFormulation
Various Methods and Formulations available.
Definition: Formulas.h:26
@ NCAR
NCAR specific formulation.
Definition: Formulas.h:29
@ DEFAULT
DEFAULT formulation.
Definition: Formulas.h:31
@ LandoltBornstein
Definition: Formulas.h:36
@ NOAA
NOAA specific formulation.
Definition: Formulas.h:30
@ UKMO
UKMO specific formulation.
Definition: Formulas.h:28
int RenumberScanPosition(int scanpos)
Get renumbered scan position 1,2,3,...
Definition: Formulas.cc:372
float Height_To_Pressure_ICAO_atmos(float height, MethodFormulation formulation)
Converts height to pressure using the International Civil Aviation Organization (ICAO) atmosphere.
Definition: Formulas.cc:224
float GetWind_V(float windSpeed, float windFromDirection)
Get northward (v) wind component from wind speed and direction.
Definition: Formulas.cc:359
float GetWindSpeed(float u, float v)
Converts u and v wind component into wind speed.
Definition: Formulas.cc:337
static const float LandoltBornstein_lookuptable[]
Definition: LookupTable.h:17
real(kind_real), parameter, public deg2rad
real(kind_real), parameter, public t0c
Definition: RunCRTM.h:27
util::Duration abs(const util::Duration &duration)
static constexpr double icao_pressure_u
Definition: Constants.h:90
static constexpr double icao_lapse_rate_u
Definition: Constants.h:77
static constexpr double icao_pressure_l
Definition: Constants.h:89
static constexpr double icao_temp_surface
Definition: Constants.h:85
static constexpr double icao_height_u
Definition: Constants.h:83
static constexpr double icao_height_l
Definition: Constants.h:81
static constexpr double deg2rad
Definition: Constants.h:21
static constexpr double epsilon
Definition: Constants.h:67
static constexpr double icao_pressure_surface
Definition: Constants.h:88
static constexpr double icao_temp_isothermal_layer
Definition: Constants.h:86
static constexpr double rad2deg
Definition: Constants.h:22
static constexpr double t0c
Definition: Constants.h:24
static constexpr double mean_earth_rad
Definition: Constants.h:40
static constexpr double icao_lapse_rate_l
Definition: Constants.h:73
static constexpr double g_over_rd
Definition: Constants.h:39