WSL/SLF GitLab Repository

ProcUnventilatedT.cc 4.16 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  Copyright 2012 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
18
#include <meteoio/meteofilters/ProcUnventilatedT.h>
19
#include <cmath>
20
21
22
23
24

using namespace std;

namespace mio {

25
26
const double ProcUnventilatedT::dflt_albedo = .23;
const double ProcUnventilatedT::soil_albedo = .23; //grass
27
const double ProcUnventilatedT::snow_albedo = .85; //snow
28
29
const double ProcUnventilatedT::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo
const double ProcUnventilatedT::vw_thresh = 0.1; //wind speed threshold
30

31
32
33
34
ProcUnventilatedT::ProcUnventilatedT(const std::vector<std::string>& vec_args)
                  : ProcessingBlock("UNVENTILATED_T"), usr_albedo(dflt_albedo),
                    is_soft(true), nakamura(true)
{
35
	parse_args(vec_args);
36
	properties.stage = ProcessingProperties::first; //for the rest: default values
37
38
}

39
void ProcUnventilatedT::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
40
41
                        std::vector<MeteoData>& ovec)
{
42
	if(param!=MeteoData::TA) {
43
		stringstream ss;
44
		ss << "Can not use " << getName() << " processing on " << MeteoData::getParameterName(param);
45
46
		throw InvalidArgumentException(ss.str(), AT);
	}
47
	ovec = ivec;
48

49
	for (size_t ii=0; ii<ovec.size(); ii++){
50
		double& tmp = ovec[ii](param);
51
52
		if(tmp == IOUtils::nodata) continue; //preserve nodata values

53
54
55
		double albedo = usr_albedo;
		double iswr = ivec[ii](MeteoData::ISWR);
		const double rswr = ivec[ii](MeteoData::RSWR);
56
		const double ta = ivec[ii](MeteoData::TA);
57
58
59
60
61
62
63
64
65
66
67
68
		double vw = ivec[ii](MeteoData::VW);
		double hs = ivec[ii](MeteoData::HS);

		if(is_soft && iswr!=IOUtils::nodata && rswr!=IOUtils::nodata && rswr>5. && iswr>5.) {
			albedo = iswr / rswr;
			hs = IOUtils::nodata; //to make sure we would not try to recompute a pseudo albedo later
		}

		if(is_soft && hs!=IOUtils::nodata && iswr==IOUtils::nodata && rswr!=IOUtils::nodata) {
			if(hs>snow_thresh) iswr = snow_albedo*rswr;
			else iswr = soil_albedo*rswr;
		}
69
70
71
72
73
74
75
76
77

		if(iswr==IOUtils::nodata || ta==IOUtils::nodata || vw==IOUtils::nodata)
			continue;

		if(is_soft && 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;
		}

78
		if(vw<vw_thresh) vw = vw_thresh; //this should be around the minimum measurable wind speed on regular instruments
79
		const double rho = 1.2; // in kg/m3
80
		const double Cp = 1004.;
81
		const double X = iswr / (rho*Cp*ta*vw);
82
		if(X<1e-4) continue; //the correction does not work well for small X values
83
84
85
86
87
88
89
90
91
		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 RE = C0 + C1*X;
			tmp -= RE; //substracting the radiative error
		} else {
			const double RE = 3.1 * sqrt(X);
			tmp -= RE; //substracting the radiative error
		}
92
93
94
	}
}

95
void ProcUnventilatedT::parse_args(std::vector<std::string> vec_args) {
96
97
98
	vector<double> filter_args;

	is_soft = false;
99
	if (!vec_args.empty()){
100
		is_soft = ProcessingBlock::is_soft(vec_args);
101
102
103
104
105
106
107
	}

	convert_args(0, 1, vec_args, filter_args);

	if (filter_args.size() > 1)
		throw InvalidArgumentException("Wrong number of arguments for filter " + getName(), AT);

108
	if (filter_args.empty()){
109
		usr_albedo = dflt_albedo;
110
	} else {
111
		usr_albedo = filter_args[0];
112
	}
113
114

	nakamura = false;
115
116
117
}

} //end namespace