WSL/SLF GitLab Repository

Commit b6062187 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Some doxygen comments have been improved, A few physical constants have been...

Some doxygen comments have been improved, A few physical constants have been added. The DEM extraction from pressure fields in NETCDFIO has been improved (the previous calculation was mostly wrong).
parent e4bc6e12
......@@ -682,6 +682,7 @@ unsigned short Date::getDayOfWeek(const bool& gmt) const {
* The week number range from 1 to 53 for a leap year. The first week is the week that contains
* the first Thursday of the year. Previous days are attributed to the last week of the previous
* year (See https://en.wikipedia.org/wiki/ISO_week_date).
* @param[out] ISO_year get filled with the matching ISO year (since the fist/last few days might belong to the previous/next year)
* @param gmt convert returned value to GMT? (default: false)
*/
unsigned short Date::getISOWeekNr(int &ISO_year, const bool& gmt) const
......
......@@ -52,20 +52,32 @@ double Atmosphere::blkBody_Radiation(const double& ea, const double& T) {
}
/**
* @brief Standard atmosphere pressure
* @brief Standard atmospheric pressure as a function of the altitude.
* This is based on the following formula (with h the altitude, P<sub>0</sub>
* and T<sub>0</sub> the standard sea level pressure and temperature, L the
* dry adiabatic lapse rate and R<sub>0</sub> the earth's radius):
* \f[
* P = P_0 \cdot {\left( 1 - \frac{L \cdot R_0 \cdot h}{T_0 ( R_0 + h )} \right)}^{\frac{g}{L R}}
* \f]
* @param altitude altitude above sea level (m)
* @return standard pressure (Pa)
*/
double Atmosphere::stdAirPressure(const double& altitude) {
const double p0 = Cst::std_press; // Air and standard pressure in Pa
const double lapse_rate = 0.0065; // K m-1
const double sea_level_temp = 288.15; // K
const double expo = Cst::gravity / (lapse_rate * Cst::gaz_constant_dry_air);
const double R0 = Cst::earth_R0; // Earth's radius in m
const double p = p0 * pow( 1. - ( (lapse_rate * R0 * altitude) / (sea_level_temp * (R0 + altitude)) ), expo );
const double expo = Cst::gravity / (Cst::dry_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
const double p = Cst::std_press * pow( 1. - ( (Cst::dry_adiabatique_lapse_rate * Cst::earth_R0 * altitude) / (Cst::std_temp * (Cst::earth_R0 + altitude)) ), expo );
return p;
}
return(p);
/**
* @brief Acceleration due to gravity
* @param altitude altitude above sea level (m)
* @param latitude latitude in degrees
* @return acceleration due to gravity (m/s2)
*/
double Atmosphere::gravity(const double& altitude, const double& latitude) {
const double lat = latitude*Cst::to_rad;
const double g = 9.780356 * (1. + 0.0052885*Optim::pow2(sin(lat)) - 0.0000059*Optim::pow2(sin(2.*lat))) - 0.003086 * altitude*1e-3;
return g;
}
/**
......
......@@ -32,6 +32,7 @@ namespace mio {
*/
class Atmosphere {
public:
static double gravity(const double& altitude, const double& latitude);
static double stdAirPressure(const double& altitude);
static double stdDryAirDensity(const double& altitude, const double& temperature);
static double waterSaturationPressure(const double& T);
......
......@@ -25,7 +25,8 @@ namespace mio {
namespace Cst {
const double stefan_boltzmann = 5.670373e-8; // (W m-2 K-4)
const double gravity = 9.80665; // (m s-2)
const double std_press = 101325.; // (Pa)
const double std_press = 101325.; // (Pa) at sea level
const double std_temp = 288.15; // (K) at sea level
const double dry_adiabatique_lapse_rate = 0.0065; // (K/m)
const double gaz_constant_dry_air = 287.058; // (J kg-1 K-1)
......@@ -35,11 +36,11 @@ namespace Cst {
const double p_water_triple_pt = 611.73; // (Pa)
const double t_water_freezing_pt = 273.15; // (K)
const double t_water_triple_pt = 273.16; // (K)
const double l_water_sublimation = 2.838e6; // (J Kg-1)
const double l_water_vaporization = 2.504e6; // (J Kg-1)
const double l_water_sublimation = 2.838e6; // (J kg-1)
const double l_water_vaporization = 2.504e6; // (J kg-1)
const double l_water_fusion = 3.34e5; // (J Kg-1)
const double water_molecular_mass = 18.0153e-3; // (Kg)
const double dry_air_mol_mass =0.0289644 ; // (Kg/mol)
const double water_molecular_mass = 18.0153e-3; // (kg)
const double dry_air_mol_mass =0.0289644 ; // (kg/mol)
const double specific_heat_ice = 2100.0; // (J K-1), at 0C
const double specific_heat_water = 4190.0; // (J K-1) at 0C
......
......@@ -269,10 +269,17 @@ void NetCDFIO::readDEM(DEMObject& dem_out)
read2DGrid_internal(p, filename, MeteoGrids::P, Date(2015,1,1,12,0,0)) &&
read2DGrid_internal(ta, filename, MeteoGrids::TA, Date(2015,1,1,12,0,0))) {
dem_out.set(p, IOUtils::nodata);
const double k = Cst::gravity*Cst::dry_air_mol_mass / (Cst::gaz_constant*Cst::dry_adiabatique_lapse_rate);
/*const double k = Cst::gravity*Cst::dry_air_mol_mass / (Cst::gaz_constant*Cst::dry_adiabatique_lapse_rate);
const double k_inv = 1./k;
for(size_t ii=0; ii<(dem_out.getNx()*dem_out.getNy()); ii++)
dem_out(ii) = ta(ii)/Cst::dry_adiabatique_lapse_rate * (pow(p(ii)/p_sea(ii), -k_inv) - 1.);
return;*/
const double k = Cst::gravity / (Cst::dry_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
const double k_inv = 1./k;
for(size_t ii=0; ii<(dem_out.getNx()*dem_out.getNy()); ii++) {
const double K = pow(p(ii)/p_sea(ii), k_inv);
dem_out(ii) = ta(ii)*Cst::earth_R0*(1.-K) / (Cst::dry_adiabatique_lapse_rate * Cst::earth_R0 - ta(ii)*(1.-K));
}
return;
}
}
......
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