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

Finished implementation of hypsometric equation.

parent ee2ccd1b
Pipeline #12618 failed with stages
in 4 minutes and 28 seconds
Hide whitespace changes
Inline Side-by-side
radixio/gfsfile.cc
View file @ 75912be7
......@@ -716,13 +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;
float relHum = 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)
......@@ -799,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, relHum, temp, msle) - surfaceAltitude;
}
// check for time
auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
if (timeIt != columns.end())
......@@ -818,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;
......@@ -840,10 +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], relHum, 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
......
radixmath/tests/tstUtil.cc
View file @ 75912be7
......@@ -11,7 +11,7 @@ TEST(radix, hpaToAltitude)
std::vector<double> pres;
std::vector<double> relhum;
std::vector<double> temp;
alt = {20, 39, 260, 294, 356, 483, 639,
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,
......@@ -28,7 +28,7 @@ TEST(radix, hpaToAltitude)
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,
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,
......@@ -40,7 +40,7 @@ TEST(radix, hpaToAltitude)
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,
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,
......@@ -52,32 +52,40 @@ TEST(radix, hpaToAltitude)
-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};
double true_altitude = 0;
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, relhumidity, temperature);
EXPECT_DOUBLE_EQ(blessed_altitude, altitude);
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.0);
}
}
......
radixmath/util.cc
View file @ 75912be7
......@@ -8,6 +8,8 @@
#include "radixmath/util.hh"
#include <algorithm>
#include <cstddef>
#include <iomanip>
#include "radixmath/point3d.hh"
#include "radixmath/vector3d.hh"
......@@ -156,25 +158,6 @@ Real cspanf(Real value, Real begin, Real end)
}
}
Real hpaToAltitude(Real hpa, Real relHum, Real temperature, Real msle)
{
// convert temperature to virtual temperature
temperature = virtualTemperature(temperature, hpa, relHum);
Real result = (SPECIFIC_GAS_CONSTANT * temperature) /
GRAVITATIONAL_ACCELERATION * std::log(msle / hpa);
// divide out the kg/mol using the molar mass of dry air and water
result = result / MOLAR_MASS_DRY_AIR;
radix_tagged_line(result << " = hpaToAltitude(" << hpa << ", " << temperature
<< ", " << msle << ")");
return result;
}
/** Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1)
{
}*/
double gammaRayAbsorptionInAir(double energy, bool linear)
{
// supported energies in MeV
......@@ -308,17 +291,4 @@ 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 @ 75912be7
......@@ -18,34 +18,62 @@
namespace radix
{
/**
* @brief mixingRatio Calculate the water vaper mixing ratio
* @param temperature kelvin
* @param pressure hPa
* @param relHum %
* @return
* @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
*/
Real mixingRatio(Real temperature, Real pressure, Real relHum);
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 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
* @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
*/
Real virtualTemperature(Real temperature, Real pressure, Real relHum);
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;
}
/**
* @brief hpaToAltitude converts hectopascals or millibars to altitude (meters)
* @param hpa hectorpascals in millibars
* @param relHum in %
* @param temperature in kelvin
* @param msle mean sea level pressure (default=1013.25)
* @return altitude in meters
*/
Real RADIX_PUBLIC hpaToAltitude(Real hpa, Real relHum, Real temperature,
Real msle = 1013.25);
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
......
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