Skip to content
GitLab
Commit f5efa905 authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Updating hpaToAltitude to hypsometric formula.

parent 273c4b86
Pipeline #12604 failed with stages
in 4 minutes and 32 seconds
Hide whitespace changes
Inline Side-by-side
radixio/gfsfile.cc
View file @ f5efa905
......@@ -718,6 +718,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
// SOUND section of Fortran
float tpot = 0.0f;
float temp = 0.0f;
float relHum = 0.0f;
bool sfcwnd = false;
float offset = 0.0f;
float plevel = 0.0f;
......@@ -802,8 +803,12 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
auto hIt = std::find(columns.begin(), columns.end(), "HGTS");
if (hIt != columns.end())
{
float v_temp = radix::virtualTemperature(temp, plevel, relHum);
radix_tagged_line("Converting temp(" << temp << ") to vtemp(" << v_temp
<< ")");
results[ll][hIt - columns.begin()] =
hpaToAltitude(plevel, temp, msle) - surfaceAltitude;
hpaToAltitude(plevel, v_temp, msle) - surfaceAltitude;
}
// check for time
auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
......@@ -841,7 +846,10 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
if (varb.compare("RH2M") == 0) varb = "RELH";
if (varb.compare("PRSS") == 0)
{
surfaceAltitude = hpaToAltitude(vdata[kk][ll], temp, msle) - 2.f;
float v_temp = radix::virtualTemperature(temp, vdata[kk][ll], relHum);
radix_tagged_line("Converting temp(" << temp << ") to vtemp(" << v_temp
<< ")");
surfaceAltitude = hpaToAltitude(vdata[kk][ll], v_temp, msle) - 2.f;
}
// check whether relative humidity has been requested, and convert
// specific humidity if this is present instead
......@@ -851,8 +859,9 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
varb = "RELH";
double specificHumidity = vdata[kk][ll];
// Convert specific humidity to relative humidity
vdata[kk][ll] =
relHum =
specificHumidityToRelativeHumidity(specificHumidity, temp, plevel);
vdata[kk][ll] = relHum;
}
auto it = std::find(columns.begin(), columns.end(), varb);
......
radixmath/constants.hh
View file @ f5efa905
......@@ -41,6 +41,15 @@ const Real SECONDS_PER_HOUR = 3600.0;
const Real METERS_PER_MILE = 1609.344;
const Real ABS_ZERO_CELSIUS = -273.15;
const Real MASS_RATIO_WATER_VAPOR_DRY_AIR = 0.62198;
/**
* @brief SPECIFIC_GAS_CONSTANT J K^-1 mol^-1
*/
const Real SPECIFIC_GAS_CONSTANT = 8.3144598;
/**
* @brief GRAVITATIONAL_ACCELERATION m/s^2
*/
const Real GRAVITATIONAL_ACCELERATION = 9.80616;
/**
* Type def unsigned long to Identifier
*/
......
radixmath/tests/tstUtil.cc
View file @ f5efa905
......@@ -7,10 +7,78 @@ using namespace radix;
TEST(radix, hpaToAltitude)
{
Real altitude = hpaToAltitude(835);
EXPECT_NEAR(1601.9551358843219, altitude, 1e-1);
altitude = hpaToAltitude(67.7);
EXPECT_NEAR(17835.433862090933, altitude, 1e-1);
std::vector<double> alt;
std::vector<double> pres;
std::vector<double> relhum;
std::vector<double> temp;
alt = {20, 39, 260, 294, 356, 483, 639,
684, 759, 1087, 1214, 1334, 1382, 1614,
1820, 1914, 2059, 2260, 2353, 2439, 2578,
2620, 2716, 2954, 3260, 3431, 3560, 3698,
3959, 4100, 4234, 4398, 4569, 5073, 5205,
5463, 5554, 5968, 6063, 6284, 6423, 6894,
7181, 7250, 7347, 7459, 7553, 7630, 7740,
7942, 8278, 8609, 9151, 9221, 9571, 9748,
1.02E+04, 1.03E+04, 1.06E+04, 1.07E+04, 1.08E+04, 1.09E+04, 1.10E+04,
1.15E+04, 1.17E+04, 1.18E+04, 1.21E+04, 1.23E+04, 1.23E+04, 1.25E+04,
1.28E+04, 1.30E+04, 1.34E+04, 1.36E+04, 1.37E+04, 1.39E+04, 1.40E+04,
1.43E+04, 1.47E+04, 1.51E+04, 1.53E+04, 1.58E+04, 1.61E+04, 1.63E+04,
1.65E+04, 1.70E+04, 1.85E+04, 1.85E+04, 1.96E+04, 2.05E+04, 2.07E+04,
2.09E+04, 2.13E+04, 2.24E+04, 2.26E+04, 2.35E+04, 2.40E+04, 2.41E+04,
2.54E+04, 2.56E+04, 2.66E+04, 2.71E+04, 2.75E+04, 2.82E+04, 2.89E+04,
2.94E+04, 2.98E+04, 3.07E+04, 3.11E+04, 3.12E+04, 3.15E+04, 3.18E+04,
3.26E+04, 3.30E+04, 3.36E+04};
pres = {1002, 1000, 973.4, 969.4, 962.1, 947.7, 930, 925, 916.6, 881,
867.5, 854.9, 850, 826.2, 805.8, 796.5, 782.4, 763.2, 754.5, 746.5,
733.7, 729.8, 721.1, 700, 673.3, 658.9, 648.2, 636.9, 616, 604.9,
594.5, 582.1, 569.3, 532.9, 523.6, 506, 500, 472.9, 466.8, 453,
444.4, 416.4, 400, 396.1, 390.7, 384.5, 379.4, 375.2, 369.3, 358.7,
341.6, 325.3, 300, 296.8, 281.4, 273.9, 253.8, 250, 242.1, 237.6,
232.3, 230, 226.6, 210.9, 201.6, 200, 190.2, 185.5, 183.6, 180,
171.4, 166, 156.3, 150.7, 150, 143.9, 141.1, 135.3, 126.7, 119.9,
115.6, 107.4, 102.8, 100, 96.8, 89.7, 70.6, 70, 59.4, 51.2,
50, 48.2, 45.5, 38.4, 37, 32.1, 30, 29.1, 24, 23.3,
20, 18.3, 17.2, 15.6, 14, 12.9, 12.1, 10.5, 10, 9.9,
9.4, 9, 8, 7.6, 7};
temp = {8.3, 7.9, 6.8, 7.4, 9.9, 10.1, 9.8, 9.5, 9.1, 6.9,
7.3, 7.8, 7.5, 6.5, 4.9, 4.4, 3.3, 1.9, 1.2, 0.8,
0, -0.2, -0.2, -2.4, -5.1, -4.2, -4.6, -5.8, -7.6, -8.1,
-9.1, -9.8, -10.5, -14.8, -15, -16.6, -17.4, -20.7, -21.4, -23.3,
-24.6, -28.4, -31.2, -31.8, -32.3, -33.2, -34, -34.6, -35.4, -36.9,
-40, -42.6, -47.4, -48.1, -49.4, -49.3, -53, -52.7, -53.6, -50,
-50.9, -49.8, -49.5, -52.8, -51.7, -52, -54, -52.7, -53, -50.8,
-52.5, -51.6, -53.4, -51.4, -51.5, -52.6, -51.7, -51.7, -52.7, -54.7,
-54.1, -55.6, -55.4, -54.1, -53.3, -56.4, -53.3, -53.5, -55.5, -54.9,
-56, -57.1, -54.2, -55, -53.5, -55.9, -54.5, -53.3, -55.5, -53.9,
-53.2, -53.6, -51, -50.9, -46, -47.2, -44.1, -41.5, -42.1, -42.7,
-40.5, -41, -41, -41.5, -38.5};
relhum = {
97.99, 98.65, 99.32, 99.32, 94.14, 77.27, 68.13, 64.8, 59.83,
64.2, 58.95, 67.71, 69.6, 75.07, 82.09, 76.34, 84.28, 79.33,
87.15, 96.46, 96.43, 95.72, 77.26, 83.58, 94.11, 73.11, 61.9,
67.25, 62.14, 49.19, 55.44, 40.29, 40.42, 47.11, 23.1, 19.71,
23.03, 28.22, 38.24, 37.59, 66.82, 55.77, 64.47, 64.95, 47.72,
47.41, 33.46, 33.61, 53.55, 53.63, 57.42, 40.18, 38.36, 38.09,
19.56, 12.76, 14.09, 13.22, 11.23, 8.743, 7.284, 6.51, 5.359,
4.952, 2.912, 2.659, 2.01, 1.508, 1.437, 1.356, 1.288, 1.065,
1.056, 0.9242, 0.984, 1.011, 0.8786, 0.8786, 0.9067, 1.02, 0.9654,
0.9886, 0.9819, 1.017, 0.9243, 1.034, 0.8928, 0.9466, 0.9766, 0.9405,
1.002, 1.05, 0.9272, 0.9689, 0.9143, 1.025, 0.9613, 0.9405, 1.048,
0.959, 0.9621, 0.9747, 0.9284, 0.9656, 0.8985, 0.9002, 0.9012, 0.9088,
0.9092, 0.9096, 0.9131, 0.9181, 0.9181, 0.9088, 0.922};
for (size_t i = 0; i < alt.size(); ++i)
{
double blessed_altitude = alt[i];
double temperature = temp[i] + 273.15;
double pressure = pres[i];
double relhumidity = relhum[i];
// calculate virtual temperature
temperature = virtualTemperature(temperature, pressure, relhumidity);
double altitude = hpaToAltitude(pressure, temperature);
EXPECT_DOUBLE_EQ(blessed_altitude, altitude);
}
}
TEST(radix, cspanf)
......
radixmath/util.cc
View file @ f5efa905
......@@ -158,7 +158,12 @@ Real cspanf(Real value, Real begin, Real end)
Real hpaToAltitude(Real hpa, Real temperature, Real msle)
{
return ((1 - std::pow(hpa / msle, 0.19022256039)) * temperature) / 0.0065;
// return ((1 - std::pow(hpa / msle, 0.19022256039)) * temperature) / 0.0065;
Real result = (SPECIFIC_GAS_CONSTANT * temperature) /
GRAVITATIONAL_ACCELERATION * std::log(msle / hpa);
radix_tagged_line(result << " = hpaToAltitude(" << hpa << ", " << temperature
<< ", " << msle << ")");
return result;
}
/** Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1)
......@@ -299,4 +304,17 @@ Real specificHumidityToRelativeHumidity(Real specificHumidity, Real absTemp,
return relHum;
}
Real virtualTemperature(Real temperature, Real pressure, Real relHum)
{
Real mR = mixingRatio(temperature, pressure, relHum);
return temperature * (1. + (mR / MASS_RATIO_WATER_VAPOR_DRY_AIR)) / (1. + mR);
}
Real mixingRatio(Real temperature, Real pressure, Real relHum)
{
Real e_s = saturationVaporPressure(temperature, pressure);
Real e_a = e_s * (relHum / 100.);
return (0.622 * e_a) / (pressure - e_a);
}
} // namespace radix
radixmath/util.hh
View file @ f5efa905
......@@ -17,6 +17,25 @@
namespace radix
{
/**
* @brief mixingRatio Calculate the water vaper mixing ratio
* @param temperature kelvin
* @param pressure hPa
* @param relHum %
* @return
*/
Real mixingRatio(Real temperature, Real pressure, Real relHum);
/**
* @brief virtualTemperature The temperature of dry air so that the its density
* is equivalent to a sample of moist air at the same pressure.
* @param temperature kelvin
* @param pressure hPa
* @param relHum %
* @return
*/
Real virtualTemperature(Real temperature, Real pressure, Real relHum);
/**
* @brief hpaToAltitude converts hectopascals or millibars to altitude (meters)
* @param hpa hectorpascals in millibars
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment