12 #include "eckit/exception/Exceptions.h"
13 #include "oops/util/Logger.h"
24 if (input ==
"UKMO")
return UKMO;
25 if (input ==
"NCAR")
return NCAR;
26 if (input ==
"NOAA")
return NOAA;
37 if (input ==
"Sonntag")
return Sonntag;
38 if (input ==
"Walko")
return Walko;
39 if (input ==
"Murphy")
return Murphy;
46 const float missingValueFloat = util::missingValue(1.0f);
48 float e_sub_s = missingValueFloat;
50 switch (formulation) {
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));
75 if (temp_K != missingValueFloat) {
80 const float Low_temp_thd = 183.15;
81 const float High_temp_thd = 338.15;
82 const float Delta_Temp = 0.1;
85 adj_Temp = std::max(Low_temp_thd, temp_K);
86 adj_Temp = std::min(High_temp_thd, adj_Temp);
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) *
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])))))));
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));
124 e_sub_s = 1000. * 0.6112 * std::exp(17.67f * (temp_K -
t0c) / (temp_K - 29.65f));
136 switch (formulation) {
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;
154 std::string errString =
"Aborting, no method matches enum formulas::MethodFormulation";
155 oops::Log::error() << errString;
156 throw eckit::BadValue(errString);
164 float QSat = util::missingValue(1.0f);
167 switch (formulation) {
188 float Tv = util::missingValue(1.0f);
190 switch (formulation) {
206 float Tv = util::missingValue(1.0f);
208 switch (formulation) {
213 float PsatRh = Psat * Rh * 0.01f;
225 const float missingValueFloat = util::missingValue(1.0f);
226 float Pressure = missingValueFloat;
228 switch (formulation) {
233 float RepT_Bot, RepT_Top, ZP1, ZP2;
240 if (height <= missingValueFloat) {
241 Pressure = missingValueFloat;
242 }
else if (height < -5000.0) {
245 Pressure = missingValueFloat;
249 Pressure = std::pow((1.0 - Pressure), ZP1);
255 Pressure = 100.0 * std::exp(Pressure);
261 std::pow((1.0 - Pressure), ZP2);
272 const float missingValueFloat = util::missingValue(1.0f);
273 const float pressure_hPa = pressure * 0.01;
274 float height = missingValueFloat;
280 if (pressure == missingValueFloat || pressure <= 0.0f) {
281 height = missingValueFloat;
282 }
else if (pressure_hPa <= 120.0f &&
293 height = 44307.692 * (1.0 - std::pow(pressure / 101325.0, 0.190));
300 if (pressure == missingValueFloat || pressure <= 0.0f) {
301 height = missingValueFloat;
324 const float missing = util::missingValue(1.0f);
328 if (u != 0.0f || v != 0.0f) {
331 windDirection = 0.0f;
334 return windDirection;
338 const float missing = util::missingValue(1.0f);
342 windSpeed = hypot(u, v);
347 float GetWind_U(
float windSpeed,
float windFromDirection) {
348 const float missing = util::missingValue(1.0f);
351 if (windFromDirection !=
missing
352 && windSpeed !=
missing && windSpeed >= 0) {
359 float GetWind_V(
float windSpeed,
float windFromDirection) {
360 const float missing = util::missingValue(1.0f);
363 if (windFromDirection !=
missing
364 && windSpeed !=
missing && windSpeed >= 0) {
374 int newpos = (scanpos + 1)/3;
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,
393 const float missingValueFloat = util::missingValue(1.0f);
395 switch (formulation) {
401 const size_t loc0 = locs.front();
404 const double lat0 = lat_in[loc0];
405 const double lon0 = lon_in[loc0];
406 const util::DateTime time0 = time_in[loc0];
412 std::vector<size_t> locs_valid;
413 for (
size_t jloc : locs) {
415 if (!apply[jloc])
continue;
417 if (windspd[jloc] != missingValueFloat && height[jloc] != missingValueFloat)
418 locs_valid.push_back(jloc);
422 if (locs_valid.size() < 2)
return;
425 const double ascent_speed = 5.16;
431 double dt_cumul = 0.0;
433 for (
size_t k = 0; k < locs_valid.size() - 1; ++k) {
435 const size_t loc_current = locs_valid[k];
436 const size_t loc_next = locs_valid[k + 1];
440 const double dh = height[loc_next] - height[loc_current];
442 const double dt = dh / ascent_speed;
446 const double avgu = 0.5 *
447 (windspd[loc_current] * std::sin((winddir[loc_current] - 180.0) *
Constants::deg2rad) +
449 const double avgv = 0.5 *
450 (windspd[loc_current] * std::cos((winddir[loc_current] - 180.0) *
Constants::deg2rad) +
461 lat_out[loc_next] = lat_out[loc_current] + dlat;
462 lon_out[loc_next] = lon_out[loc_current] + dlon;
465 time_out[loc_next] = time0 + util::Duration(
static_cast<int64_t
>(dt_cumul));
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()) {
477 time = time_out[jloc];
481 time_out[jloc] = time;
static const float LandoltBornstein_lookuptable[]
real(kind_real), parameter, public deg2rad
real(kind_real), parameter, public t0c
util::Duration abs(const util::Duration &duration)
static constexpr double icao_pressure_u
static constexpr double icao_lapse_rate_u
static constexpr double icao_pressure_l
static constexpr double icao_temp_surface
static constexpr double icao_height_u
static constexpr double icao_height_l
static constexpr double deg2rad
static constexpr double epsilon
static constexpr double icao_pressure_surface
static constexpr double icao_temp_isothermal_layer
static constexpr double rad2deg
static constexpr double t0c
static constexpr double mean_earth_rad
static constexpr double icao_lapse_rate_l
static constexpr double g_over_rd