WSL/SLF GitLab Repository

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

Full rewrite of the Huwald and Nakamura algorithms. Now they perform properly!

parent c893b9cb
......@@ -56,46 +56,60 @@ void ProcUnventilatedT::filterTA(const unsigned int& param, std::vector<MeteoDat
void ProcUnventilatedT::correctTA(const unsigned int& param, std::vector<MeteoData>& ovec) const
{
for (size_t ii=0; ii<ovec.size(); ii++) {
double& ta = ovec[ii](param);
if (ta == IOUtils::nodata) continue; //preserve nodata values
double albedo = usr_albedo;
double iswr = ovec[ii](MeteoData::ISWR);
const double rswr = ovec[ii](MeteoData::RSWR);
double vw = ovec[ii](MeteoData::VW);
double hs = ovec[ii](MeteoData::HS);
if (iswr!=IOUtils::nodata && rswr!=IOUtils::nodata && rswr>5. && iswr>5.) {
albedo = rswr / iswr;
hs = IOUtils::nodata; //to make sure we would not try to recompute a pseudo albedo later
}
if (hs!=IOUtils::nodata && iswr==IOUtils::nodata && rswr!=IOUtils::nodata) {
if (hs>snow_thresh) iswr = snow_albedo*rswr;
else iswr = soil_albedo*rswr;
}
if (iswr==IOUtils::nodata || vw==IOUtils::nodata)
continue;
if (nakamura) {
for (size_t ii=0; ii<ovec.size(); ii++) {
double& ta = ovec[ii](param);
if (ta == IOUtils::nodata) continue; //preserve nodata values
double iswr = ovec[ii](MeteoData::ISWR);
const double rswr = ovec[ii](MeteoData::RSWR);
double vw = ovec[ii](MeteoData::VW);
const double hs = ovec[ii](MeteoData::HS);
if (hs!=IOUtils::nodata && iswr==IOUtils::nodata && rswr!=IOUtils::nodata) {
if (hs>snow_thresh) iswr = rswr / snow_albedo;
else iswr = rswr / soil_albedo;
}
if (iswr==IOUtils::nodata || vw==IOUtils::nodata)
continue;
if (vw<vw_thresh) vw = vw_thresh; //this should be around the minimum measurable wind speed on regular instruments
const double rho = 1.2; // in kg/m3
const double Cp = 1004.;
const double X = iswr / (rho*Cp*ta*vw);
if (X<1e-4) continue; //the correction does not work well for small X values
if (hs!=IOUtils::nodata) { //try to get snow height in order to adjust the albedo
if (hs>snow_thresh) albedo = snow_albedo;
else albedo = soil_albedo;
}
if (vw<vw_thresh) vw = vw_thresh; //this should be around the minimum measurable wind speed on regular instruments
const double rho = 1.2; // in kg/m3
const double Cp = 1004.;
const double X = iswr / (rho*Cp*ta*vw);
if (X<1e-4) continue; //the correction does not work well for small X values
if (nakamura) {
const double C0 = 0.13;
const double C1 = 373.40; /* albedo / dflt_albedo; //in order to introduce the albedo as a scaling factor*/
const double C1 = 373.40;
const double RE = C0 + C1*X;
ta -= RE; //substracting the radiative error
} else {
const double RE = 3.1 * sqrt(albedo*X*1e3);
}
} else { //Huwald
for (size_t ii=0; ii<ovec.size(); ii++) {
double& ta = ovec[ii](param);
if (ta == IOUtils::nodata) continue; //preserve nodata values
const double iswr = ovec[ii](MeteoData::ISWR);
double rswr = ovec[ii](MeteoData::RSWR);
double vw = ovec[ii](MeteoData::VW);
const double hs = ovec[ii](MeteoData::HS);
if (hs!=IOUtils::nodata && rswr==IOUtils::nodata && iswr!=IOUtils::nodata) {
if (hs>snow_thresh) rswr = iswr * snow_albedo;
else rswr = iswr * soil_albedo;
}
if (rswr==IOUtils::nodata || vw==IOUtils::nodata)
continue;
if (vw<vw_thresh) vw = vw_thresh; //this should be around the minimum measurable wind speed on regular instruments
const double rho = 1.2; // in kg/m3
const double Cp = 1004.;
const double X = rswr / (rho*Cp*ta*vw) * 1e3;
if (X<0.05) continue; //the correction does not work well for small X values
const double RE = 3.1 * sqrt(X);
ta -= RE; //substracting the radiative error
}
}
......
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