WSL/SLF GitLab Repository

ProcUndercatch_WMO.cc 5.99 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/ProcUndercatch_WMO.h>
19
#include <meteoio/meteolaws/Atmosphere.h>
20
21
22
23
24
25
#include <cmath>

using namespace std;

namespace mio {

26
const double ProcUndercatch_WMO::Tsnow_WMO=-2., ProcUndercatch_WMO::Train_WMO=2.; //WMO values from Yan et al (2001)
27

28
ProcUndercatch_WMO::ProcUndercatch_WMO(const std::vector<std::string>& vec_args) : ProcessingBlock("UNDERCATCH_WMO") {
29
30
	Tsnow = Tsnow_WMO;
	Train = Train_WMO;
31
32
33
34
	parse_args(vec_args);
	properties.stage = ProcessingProperties::first; //for the rest: default values
}

35
void ProcUndercatch_WMO::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
36
37
                        std::vector<MeteoData>& ovec)
{
38
	if(param!=MeteoData::HNW)
39
		throw InvalidArgumentException("Trying to use UNDERCATCH_WMO filter on " + MeteoData::getParameterName(param) + " but it can only be applied to precipitation!!" + getName(), AT);
40
41
42
	ovec.clear();
	ovec.reserve(ivec.size());

43
	for (size_t ii=0; ii<ivec.size(); ii++){
44
45
		ovec.push_back(ivec[ii]);

46
		double& tmp = ovec[ii](param);
47
		double VW = ovec[ii](MeteoData::VW);
48
		if(VW!=IOUtils::nodata) VW = Atmosphere::windLogProfile(VW, 10., 2.); //impact seems minimal
49
		double t = ovec[ii](MeteoData::TA);
50
51
		if(t==IOUtils::nodata) continue; //we MUST have air temperature in order to filter
		t=K_TO_C(t); //t in celsius
52
53
54
55
56
57
58
59
60
61
62
63
		precip_type precip=mixed;
		if(t<=Tsnow) precip=snow;
		if(t>=Train) precip=rain;

		//HACK:: we don't use Tmax, Tmin, Tmean but only the current temperature instead
		if (tmp == IOUtils::nodata || tmp==0. || precip==rain) {
			continue; //preserve nodata values and no precip or purely liquid precip
		} else if(type==cst) {
			if(precip==snow) tmp *= factor_snow;
			if(precip==mixed) tmp *= factor_mixed;
		} else if(type==nipher) {
			if(VW==IOUtils::nodata) continue;
64
			double k=100.;
65
66
67
68
69
70
71
72
			if(precip==snow) k=100.-0.44*VW*VW-1.98*VW;
			if(precip==mixed) {
				if(t==IOUtils::nodata) continue;
				k=97.29-3.18*VW+0.58*t-0.67*t; //Tmax, Tmin
			}
			tmp *= 100./k;
		} else if(type==tretyakov) {
			if(VW==IOUtils::nodata || t==IOUtils::nodata) continue;
73
			double k=100.;
74
			if(VW>8.5) VW=8.5; //the fits have been calibrated until 8.5 m/s
75
76
77
78
79
			if(precip==snow) k=103.11-8.67*VW+0.30*t; //Tmax
			if(precip==mixed) k=96.99-4.46*VW+0.88*t+0.22*t; //Tmax, Tmin
			tmp *= 100./k;
		} else if(type==us8sh) {
			if(VW==IOUtils::nodata) continue;
80
			double k=100.;
81
82
83
84
85
			if(precip==snow) k=exp(4.61-0.04*pow(VW, 1.75));
			if(precip==mixed) k=101.04-5.62*VW;
			tmp *= 100./k;
		} else if(type==us8unsh) {
			if(VW==IOUtils::nodata) continue;
86
			double k=100.;
87
88
89
			if(precip==snow) k=exp(4.61-0.16*pow(VW, 1.28));
			if(precip==mixed) k=100.77-8.34*VW;
			tmp *= 100./k;
90
		} else if(type==rt3_jp) {
91
92
93
			if(VW==IOUtils::nodata) continue;
			const double rh = ovec[ii](MeteoData::RH);
			const double alt = ovec[ii].meta.position.getAltitude();
94
			double k=100.;
95
			if(rh!=IOUtils::nodata && alt!=IOUtils::nodata) {
96
				const double t_wb = K_TO_C(Atmosphere::wetBulbTemperature(ovec[ii](MeteoData::TA), rh, alt));
97
98
99
100
101
				double ts_rate;
				if(t_wb<1.1) ts_rate = 1. - .5*exp(-2.2*pow(1.1-t_wb, 1.3));
				else ts_rate = .5*exp(-2.2*pow(t_wb-1.1, 1.3));
				if(ts_rate>.5) precip=snow; else precip=mixed;
			}
102
103
			if(precip==snow) k=100. / (1.+.346*VW);
			if(precip==mixed) k=100. / (1.+.0856*VW);
104
			tmp *= 100./k;
105
106
		} else if(type==hellmann) {
			if(VW==IOUtils::nodata) continue;
107
			double k=100.;
108
109
110
111
112
113
			if(precip==snow) k=100.+1.13*VW*VW-19.45*VW;
			if(precip==mixed) {
				if(t==IOUtils::nodata) continue;
				k=96.63+0.41*VW*VW-9.84*VW+5.95*t; //Tmean
			}
			tmp *= 100./k;
114
115
		}  else if(type==hellmannsh) {
			if(VW==IOUtils::nodata) continue;
116
			double k=100.;
117
118
119
			if(precip==snow) k=100.+0.72*VW*VW-13.74*VW;
			if(precip==mixed) k=101.319+0.524*VW*VW-6.42*VW;
			tmp *= 100./k;
120
121
122
123
		}
	}
}

124
void ProcUndercatch_WMO::parse_args(std::vector<std::string> filter_args)
125
126
127
128
129
130
131
132
133
134
{
	if (filter_args.size() < 1)
		throw InvalidArgumentException("Wrong number of arguments for filter " + getName(), AT);

	for(size_t ii=0; ii<filter_args.size(); ii++) {
		IOUtils::toLower(filter_args[ii]);
	}

	if(filter_args[0]=="cst") {
		type=cst;
135
		if(filter_args.size() < 3 || filter_args.size() > 5 || filter_args.size() == 4)
136
137
138
			throw InvalidArgumentException("Wrong number of arguments for filter "+getName()+" with rain gauge type \"cst\"", AT);
		IOUtils::convertString(factor_snow, filter_args[1]);
		IOUtils::convertString(factor_mixed, filter_args[2]);
139
140
141
142
143
144
		if(filter_args.size()==5) {
			IOUtils::convertString(Tsnow, filter_args[1]);
			Tsnow = K_TO_C(Tsnow);
			IOUtils::convertString(Train, filter_args[2]);
			Train = K_TO_C(Train);
		}
145
146
147
148
149
150
151
152
	} else if(filter_args[0]=="nipher") {
		type=nipher;
	} else if(filter_args[0]=="tretyakov") {
		type=tretyakov;
	} else if(filter_args[0]=="us8sh") {
		type=us8sh;
	} else if(filter_args[0]=="us8unsh") {
		type=us8unsh;
153
154
	} else if(filter_args[0]=="rt3_jp") {
		type=rt3_jp;
155
156
	} else if(filter_args[0]=="hellmann") {
		type=hellmann;
157
158
	} else if(filter_args[0]=="hellmannsh") {
		type=hellmannsh;
159
160
161
162
163
164
	} else {
		throw InvalidArgumentException("Rain gauge type \""+ filter_args[0] +"\" unknown for filter "+getName(), AT);
	}
}

} //end namespace