Skip to content

GitLab

Commit f05ea3e9 authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Merge branch 'pressure_altitude' into 'master' 

Pressure altitude

See merge request !26
parents 273c4b86 819472de
Pipeline #12623 passed with stages
in 5 minutes and 46 seconds
Hide whitespace changes
Inline Side-by-side
radixio/gfsfile.cc
View file @ f05ea3e9
......@@ -716,12 +716,13 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
delete rstream;
// SOUND section of Fortran
float tpot = 0.0f;
float temp = 0.0f;
bool sfcwnd = false;
float offset = 0.0f;
float plevel = 0.0f;
float surfaceAltitude = 0.f;
float tpot = 0.0f;
float temp = 288.15f; // default to surface temperature
float relHum = 0.0f;
bool sfcwnd = false;
float offset = 0.0f;
float plevel = 0.0f;
double surfaceAltitude = 0;
std::vector<std::vector<float>> results(mHeader.nz);
for (int ll = 0; ll < mHeader.nz; ++ll)
......@@ -798,13 +799,6 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
// initialize space for results vector
results[ll] = std::vector<float>(columns.size(), 0.0f);
// if "HGTS" has been requested add level as the default HGTS
auto hIt = std::find(columns.begin(), columns.end(), "HGTS");
if (hIt != columns.end())
{
results[ll][hIt - columns.begin()] =
hpaToAltitude(plevel, temp, msle) - surfaceAltitude;
}
// check for time
auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
if (timeIt != columns.end())
......@@ -817,7 +811,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
results[ll][timeIt - columns.begin()] = (float)hour;
}
// check for pressure
auto presIt = std::find(columns.begin(), columns.end(), "PRSS");
auto presIt = std::find(columns.begin(), columns.end(), "PRES");
if (presIt != columns.end())
{
results[ll][presIt - columns.begin()] = plevel;
......@@ -839,9 +833,12 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
// humidity(RH2M)
if (varb.compare("T02M") == 0) varb = "TEMP";
if (varb.compare("RH2M") == 0) varb = "RELH";
// convert surface specific humidity to specific humidity for relhum
// conversion
if (varb.compare("SPH2") == 0) varb = "SPHU";
if (varb.compare("PRSS") == 0)
{
surfaceAltitude = hpaToAltitude(vdata[kk][ll], temp, msle) - 2.f;
surfaceAltitude = hypsometric(msle, vdata[kk][ll], temp);
}
// check whether relative humidity has been requested, and convert
// specific humidity if this is present instead
......@@ -851,8 +848,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);
......@@ -861,18 +859,6 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
results[ll][it - columns.begin()] = vdata[kk][ll];
}
}
// check surface data (2 meters)
if (ll == 0)
{
// check for surface height
{
auto it = std::find(columns.begin(), columns.end(), "HGTS");
if (it != columns.end())
{
results[ll][it - columns.begin()] = 2.0f;
}
}
}
bool hwind = false; // have wind?
float wd = 0.f;
float ws = 0.f;
......
radixmath/constants.hh
View file @ f05ea3e9
......@@ -41,6 +41,24 @@ 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;
/**
* @brief MOLAR_MASS_DRY_AIR kg/mol
*/
const Real MOLAR_MASS_DRY_AIR = 0.028964;
/**
* @brief MOLAR_MASS_WATER_VAPER kg/mol
*/
const Real MOLAR_MASS_WATER_VAPER = 0.018016;
/**
* Type def unsigned long to Identifier
*/
......
radixmath/tests/tstUtil.cc
View file @ f05ea3e9
......@@ -7,10 +7,86 @@ 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};
double true_altitude = 0;
for (size_t i = 0; i < alt.size(); ++i)
{
double blessed_altitude = alt[i];
double temperature = temp[0] + 273.15; // surface
double prev_pres = pres[0]; // mean sea level pressure
if (0 != i)
{ // mean temperature
temperature = 273.15 + (temp[i] + temp[i - 1]) / 2.0;
prev_pres = pres[i - 1];
}
// calculates thinkness in meters between prev_pres and pres[i]
true_altitude += hypsometric(prev_pres, pres[i], temperature);
// we augment by assumed surface elevation (alt[0])
radix_tagged_line(std::setw(10) << (alt[0] + true_altitude)
<< " = hpaToAltitude(" << pres[i] << ", "
<< temperature << ", " << prev_pres << ")");
EXPECT_NEAR(blessed_altitude, (alt[0] + true_altitude), 80.2);
}
// convert temperature to kelvin
std::transform(temp.begin(), temp.end(), temp.begin(),
[](Real t) { return t + 273.15; });
// Use the convenience helper for the entire profile
std::vector<Real> hypProfile = hypsometricProfile(pres, temp, true);
// adjust hyp profile with surface elevation
std::transform(hypProfile.begin(), hypProfile.end(), hypProfile.begin(),
[&alt](Real a) { return a + alt[0]; });
// check results
for (size_t i = 0; i < alt.size(); ++i)
{
EXPECT_NEAR(alt[i], hypProfile[i], 80.2);
}
}
TEST(radix, cspanf)
......
radixmath/util.cc
View file @ f05ea3e9
......@@ -8,6 +8,8 @@
#include "radixmath/util.hh"
#include <algorithm>
#include <cstddef>
#include <iomanip>
#include "radixmath/point3d.hh"
#include "radixmath/vector3d.hh"
......@@ -156,16 +158,6 @@ 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;
}
/** Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1)
{
}*/
double gammaRayAbsorptionInAir(double energy, bool linear)
{
// supported energies in MeV
......
radixmath/util.hh
View file @ f05ea3e9
......@@ -18,14 +18,62 @@
namespace radix
{
/**
* @brief hpaToAltitude converts hectopascals or millibars to altitude (meters)
* @param hpa hectorpascals in millibars
* @param temperature in kelvin (default=288.15)
* @param msle mean sea level pressure (default=1013.25)
* @return altitude in meters
*/
Real RADIX_PUBLIC hpaToAltitude(Real hpa, Real temperature = 288.15,
Real msle = 1013.25);
* @brief hypsometric Calculates the thickness in meters between p1 and p2.
* @param p1 z1 pressure kPa
* @param p2 z2 pressure kPa
* @param meanTemp mean temperature between z1 and z2.
* @return thickness in meters
*/
template <typename type>
type RADIX_PUBLIC hypsometric(type p1, type p2, type meanTemp)
{
type result = (SPECIFIC_GAS_CONSTANT * meanTemp) /
GRAVITATIONAL_ACCELERATION * std::log(p1 / p2);
// divide out the kg/mol using the molar mass of dry air and water
result = result / MOLAR_MASS_DRY_AIR;
return result;
}
/**
* @brief hypsometricProfile Calculate thickness' for an entire profile
* @param pressure std::vector<Real> pressures in kPa or millibars
* @param temperature std::vector<Real>temperatures in Kelvin
* @param sum whether to sum each thickness so it is the thickness between 0 &
* the level.
* @return std::vector<Real> thickness in meters
*/
template <typename type>
std::vector<type> RADIX_PUBLIC
hypsometricProfile(const std::vector<type> &pressure,
const std::vector<type> &temperature, bool sum = true)
{
if (pressure.size() == 0) return std::vector<type>();
if (pressure.size() != temperature.size())
throw std::runtime_error(
"hypsometricProfile requires pressure and temperature be the same "
"size.");
std::vector<type> altitude(pressure.size(), 0);
type p1 = pressure[0];
type t1 = temperature[0];
for (size_t pi = 0; pi < pressure.size(); ++pi)
{
if (pi != 0)
{
p1 = pressure[pi - 1];
t1 = (temperature[pi] + temperature[pi - 1]) / 2.0;
}
altitude[pi] = hypsometric(p1, pressure[pi], t1);
}
if (sum)
{
for (size_t ai = 1; ai < altitude.size(); ++ai)
{
altitude[ai] = altitude[ai] + altitude[ai - 1];
}
}
return altitude;
}
/**
* @brief absoluteTemperature Get the temperature of an air parcel given its
......
Markdown is supported
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