diff --git a/Source/meteoio/CMakeLists.txt b/Source/meteoio/CMakeLists.txt
index cdc29fcc..09b8cf28 100644
--- a/Source/meteoio/CMakeLists.txt
+++ b/Source/meteoio/CMakeLists.txt
@@ -90,8 +90,8 @@ SET(PLUGIN_ALPUG ON CACHE BOOL "Compilation ALPUG ON or OFF")
 SET(PLUGIN_ARCIO ON CACHE BOOL "Compilation ARCIO ON or OFF")
 SET(PLUGIN_ARGOSIO OFF CACHE BOOL "Compilation ARGOSIO ON or OFF")
 SET(PLUGIN_ARPSIO ON CACHE BOOL "Compilation ARPSIO ON or OFF")
-SET(PLUGIN_CSVIO ON CACHE BOOL "Compilation CsvIO ON or OFF")
 SET(PLUGIN_COSMOXMLIO OFF CACHE BOOL "Compilation COSMOXMLIO ON or OFF")
+SET(PLUGIN_CSVIO ON CACHE BOOL "Compilation CsvIO ON or OFF")
 SET(PLUGIN_DBO OFF CACHE BOOL "Compilation DBO ON or OFF")
 SET(PLUGIN_GEOTOPIO ON CACHE BOOL "Compilation GEOTOPIO ON or OFF")
 SET(PLUGIN_GOESIO ON CACHE BOOL "Compilation GOESIO ON or OFF")
@@ -105,9 +105,10 @@ SET(PLUGIN_PGMIO ON CACHE BOOL "Compilation PGMIO ON or OFF")
 SET(PLUGIN_PMODIO OFF CACHE BOOL "Compilation PMODIO ON or OFF")
 SET(PLUGIN_PNGIO OFF CACHE BOOL "Compilation PNGIO ON or OFF")
 SET(PLUGIN_PSQLIO OFF CACHE BOOL "Compilation PSQLIO ON or OFF")
+SET(PLUGIN_SASEIO OFF CACHE BOOL "Compilation SASEIO ON or OFF")
 SET(PLUGIN_SMETIO ON CACHE BOOL "Compilation SMETIO ON or OFF")
 SET(PLUGIN_SNIO ON CACHE BOOL "Compilation SNIO ON or OFF")
-SET(PLUGIN_SASEIO OFF CACHE BOOL "Compilation SASEIO ON or OFF")
+SET(PLUGIN_WWCSIO OFF CACHE BOOL "Compilation WWCSIO ON or OFF")
 SET(PLUGIN_ZRXPIO OFF CACHE BOOL "Compilation ZRXPIO ON or OFF")
 SET(PROJ OFF CACHE BOOL "Use PROJ for the class MapProj ON or OFF")
 
diff --git a/Source/meteoio/doc/examples/input/surface-grids/Switzerland_1000m.asc b/Source/meteoio/doc/examples/input/surface-grids/Switzerland_1000m.asc
old mode 100755
new mode 100644
diff --git a/Source/meteoio/meteoio/Config.cc b/Source/meteoio/meteoio/Config.cc
index 3a023081..6f4256bd 100644
--- a/Source/meteoio/meteoio/Config.cc
+++ b/Source/meteoio/meteoio/Config.cc
@@ -269,11 +269,12 @@ std::vector<std::string> Config::getKeys(std::string keymatch,
 	} else {
 		keymatch = section + keymatch;
 
-		for (std::map<string,string>::const_iterator it=properties.begin(); it != properties.end(); ++it) {
-			const size_t found_pos = (it->first).find(keymatch, 0);
+		//for (std::map<string,string>::const_iterator it=properties.begin(); it != properties.end(); ++it) {
+		for (const auto& it : properties) {
+			const size_t found_pos = (it.first).find(keymatch, 0);
 			if (found_pos==0) { //found it!
 				const size_t section_len = section.length();
-				const std::string key( (it->first).substr(section_len) ); //from pos to the end
+				const std::string key( (it.first).substr(section_len) ); //from pos to the end
 				vecResult.push_back( key );
 			}
 		}
@@ -291,8 +292,8 @@ void Config::write(const std::string& filename) const
 	try {
 		std::string current_section;
 		unsigned int sectioncount = 0;
-		for (std::map<string,string>::const_iterator it=properties.begin(); it != properties.end(); ++it) {
-			const std::string key_full( it->first );
+		for (const auto& it : properties) {
+			const std::string key_full( it.first );
 			const std::string section( ConfigParser::extract_section(key_full) );
 
 			if (current_section != section) {
@@ -304,7 +305,7 @@ void Config::write(const std::string& filename) const
 			}
 
 			const size_t key_start = key_full.find_first_of(":");
-			const std::string value( it->second );
+			const std::string value( it.second );
 			if (value.empty()) continue;
 
 			if (key_start!=string::npos) //start after the "::" marking the section prefix
@@ -327,7 +328,7 @@ const std::string Config::toString() const {
 	std::ostringstream os;
 	os << "<Config>\n";
 	os << "Source: " << sourcename << "\n";
-	for (map<string,string>::const_iterator it = properties.begin(); it != properties.end(); ++it){
+	for (auto it = properties.begin(); it != properties.end(); ++it){
 		os << (*it).first << " -> " << (*it).second << "\n";
 	}
 	os << "</Config>\n";
diff --git a/Source/meteoio/meteoio/Timer.cc b/Source/meteoio/meteoio/Timer.cc
index e21e92ca..f4999b1f 100644
--- a/Source/meteoio/meteoio/Timer.cc
+++ b/Source/meteoio/meteoio/Timer.cc
@@ -23,7 +23,7 @@
 	#include <windows.h>
 	#undef max
 	#undef min
-	#include <cstdlib>
+	#include <threadpoollegacyapiset.h>
 #else
 	#include <sys/time.h>
 	#include <unistd.h>
diff --git a/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.cc b/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.cc
index c108cfc9..e95ed287 100644
--- a/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.cc
+++ b/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.cc
@@ -19,11 +19,16 @@
 #include <cmath>
 #include <limits.h>
 #include <algorithm>
+#include <fstream>
+#include <sstream>
+//#include <cerrno>
+#include <cstring>
 
 #include <meteoio/dataClasses/DEMAlgorithms.h>
 #include <meteoio/dataClasses/DEMObject.h>
 #include <meteoio/MathOptim.h>
 #include <meteoio/IOUtils.h>
+#include <meteoio/FileUtils.h>
 #include <meteoio/meteoLaws/Meteoconst.h> //for math constants
 
 /**
@@ -90,7 +95,7 @@ double DEMAlgorithms::getSearchDistance(const DEMObject& dem)
 * @param[in] ix1 x index of the origin point
 * @param[in] iy1 y index of the origin point
 * @param[in] bearing direction given by a compass bearing
-* @return tangente of angle above the horizontal (in deg)
+* @return tangente of angle above the horizontal (in deg) or IOUtils::nodata if the point (ix1, iy1) does not fit within the provided DEM
 */
 double DEMAlgorithms::getHorizon(const DEMObject& dem, const size_t& ix1, const size_t& iy1, const double& bearing)
 {
@@ -100,6 +105,7 @@ double DEMAlgorithms::getHorizon(const DEMObject& dem, const size_t& ix1, const
 	
 	const int dimx = (signed)dem.grid2D.getNx();
 	const int dimy = (signed)dem.grid2D.getNy();
+	if ((signed)ix1>dimx || (signed)iy1>dimy) return IOUtils::nodata; //in case the point does not fith within the provided DEM
 	if (ix1==0 || (signed)ix1==dimx-1 || iy1==0 || (signed)iy1==dimy-1) return 0.; //a border cell is not shadded
 	
 	const double cell_alt = dem.grid2D(ix1, iy1);
@@ -135,29 +141,200 @@ double DEMAlgorithms::getHorizon(const DEMObject& dem, const size_t& ix1, const
 * @param[in] dem DEM to work with
 * @param[in] point the origin point
 * @param[in] bearing direction given by a compass bearing
-* @return tangente of angle above the horizontal (in deg)
+* @return tangente of angle above the horizontal (in deg) or IOUtils::nodata if the provided point does not fir within the 
+* provided DEM
 */
-double DEMAlgorithms::getHorizon(const DEMObject& dem, const Coords& point, const double& bearing)
+double DEMAlgorithms::getHorizon(const DEMObject& dem, Coords point, const double& bearing)
 {
-	const int ix1 = (int)point.getGridI();
-	const int iy1 = (int)point.getGridJ();
+	if (!point.indexIsValid()) {
+		if (!dem.gridify(point)) return IOUtils::nodata;
+	}
+	
+	const size_t ix1 = static_cast<size_t>( point.getGridI() );
+	const size_t iy1 = static_cast<size_t>( point.getGridJ() );
 	return getHorizon(dem, ix1, iy1, bearing);
 }
 
 /**
-* @brief Returns the horizon from a given point looking 360 degrees around by increments
+* @brief Returns the horizon from a given point looking 360 degrees around by increments. If the provided point does not
+* fit within the provided DEM, an empty result set is returned.
 * @param[in] dem DEM to work with
 * @param[in] point the origin point
 * @param[in] increment to the bearing between two angles
-* @param[out] horizon vector of heights above a given angle
-*
+* @return horizon vector of heights as a function of azimuth
 */
-void DEMAlgorithms::getHorizon(const DEMObject& dem, const Coords& point, const double& increment, std::vector< std::pair<double,double> >& horizon)
+std::vector< std::pair<double,double> > DEMAlgorithms::getHorizonScan(const DEMObject& dem, Coords point, const double& increment)
 {
+	std::vector< std::pair<double,double> > horizon;
+	if (!point.indexIsValid()) {
+		if (!dem.gridify(point)) return horizon;
+	}
+	
+	const size_t ix1 = static_cast<size_t>( point.getGridI() );
+	const size_t iy1 = static_cast<size_t>( point.getGridJ() );
 	for (double bearing=0.0; bearing <360.; bearing += increment) {
-		const double tan_alpha = getHorizon(dem, point, bearing);
-		horizon.push_back( make_pair(bearing, atan(tan_alpha)*Cst::to_deg) );
+		const double tan_alpha = getHorizon(dem, ix1, iy1, bearing);
+		if (tan_alpha!=IOUtils::nodata)
+			horizon.push_back( make_pair(bearing, atan(tan_alpha)*Cst::to_deg) );
+	}
+	
+	return horizon;
+}
+
+//custom function for sorting the horizons
+struct sort_horizons {
+	bool operator()(const std::pair<double,double> &left, const std::pair<double,double> &right) {
+		if (left.first < right.first) return true; else return false;
+	}
+};
+
+/**
+* @brief Read the horizons from a given set of points looking 360 degrees around provided in a file
+* @details The file containing the horizons is made of any number of lines with the following structure: 
+* `{stationID} {azimuth} {elevation}` where the azimuth is given in degrees North (as read from a compass) and the
+* elevation as degrees above the horizontal. For example:
+* @code
+* STB2 0 2
+* STB2 210 18
+* WFJ 180 5
+* WFJ 130 10
+* STB2 180 25
+* @endcode
+* 
+* @param[in] where Description of the caller to be used in error messages, such as 'Filter::shade'
+* @param[in] filename the file and path containing the horizon
+* @return a map containing for each stationID the horizon vector of heights as a function of azimuth 
+*/
+std::map< std::string, std::vector< std::pair<double,double> > > DEMAlgorithms::readHorizonScan(const std::string& where, const std::string& filename)
+{
+	std::ifstream fin( filename.c_str() );
+	if (fin.fail()) {
+		std::ostringstream ss;
+		ss << where << ": " << "error opening file \"" << filename << "\", possible reason: " << std::strerror(errno);
+		throw AccessException(ss.str(), AT);
+	}
+
+	const char eoln = FileUtils::getEoln(fin); //get the end of line character for the file
+	std::map< std::string, std::vector< std::pair<double,double> > > horizon;
+
+	try {
+		size_t lcount = 0;
+		double azimuth, value;
+		std::string stationID, line;
+		do {
+			lcount++;
+			getline(fin, line, eoln); //read complete line
+			IOUtils::stripComments(line);
+			IOUtils::trim(line);
+			if (line.empty()) continue;
+
+			std::istringstream iss(line);
+			iss.setf(std::ios::fixed);
+			iss.precision(std::numeric_limits<double>::digits10);
+			iss >> std::skipws >> stationID;
+			iss >> std::skipws >> azimuth;
+			if ( !iss || azimuth<0. || azimuth>360.) {
+				std::ostringstream ss;
+				ss << "Invalid azimuth in file " << filename << " at line " << lcount;
+				throw InvalidArgumentException(ss.str(), AT);
+			}
+			iss >> std::skipws >> value;
+			if ( !iss ){
+				std::ostringstream ss;
+				ss << "Invalid value in file " << filename << " at line " << lcount;
+				throw InvalidArgumentException(ss.str(), AT);
+			}
+
+			horizon[stationID].push_back( make_pair(azimuth, value) );
+		} while (!fin.eof());
+		fin.close();
+	} catch (const std::exception&){
+		if (fin.is_open()) {//close fin if open
+			fin.close();
+		}
+		throw;
+	}
+	
+	if (horizon.empty()) throw InvalidArgumentException(where+", no valid horizon found in file '"+filename+"'", AT);
+	for (auto& it : horizon) {
+		if (it.second.empty()) throw InvalidArgumentException(where+", no valid horizon for station '"+it.first+"' in file '"+filename+"'", AT);
+		std::sort(it.second.begin(), it.second.end(), sort_horizons());
+	}
+
+	return horizon;
+}
+
+/**
+* @brief Write to a file the horizons from a given set of points looking 360 degrees around
+* @details The file containing the horizons is made of any number of lines with the following structure: 
+* `{stationID} {azimuth} {elevation}` where the azimuth is given in degrees North (as read from a compass) and the
+* elevation as degrees above the horizontal. For example:
+* @code
+* STB2 0 2
+* STB2 210 18
+* WFJ 180 5
+* WFJ 130 10
+* STB2 180 25
+* @endcode
+* 
+* @param[in] horizon a map of vectors of (azimuth, elevation) coordinates defining the horizon, per stationID
+* @param[in] filename the file and path where to write the horizon
+*/
+void DEMAlgorithms::writeHorizons(const std::map< std::string, std::vector< std::pair<double,double> > >& horizon, const std::string& filename)
+{
+	if (!FileUtils::validFileAndPath(filename)) throw InvalidNameException(filename,AT);
+	std::ofstream fout(filename.c_str(), ios::out);
+	if (fout.fail()) {
+		std::ostringstream ss;
+		ss << "error opening file \"" << filename << "\" for writing, possible reason: " << std::strerror(errno);
+		throw AccessException(ss.str(), AT);
+	}
+	
+	for (const auto& it : horizon) {
+		const std::string stationID( it.first );
+		for (size_t ii=0; ii<it.second.size(); ii++) 
+			fout << stationID << " " << it.second[ii].first << " " << it.second[ii].second << "\n";
+	}
+	
+	if (fout.is_open()) //close fout if open
+		fout.close();
+}
+
+/**
+* @brief Linearly interpolate the horizon height for any given azimuth given a set of (azimuth, elevation) points
+* @param[in] horizon a set of (azimuth, elevation) coordinates defining the horizon
+* @param[in] azimuth the azimuth where to perform the interpolation
+* @return the interpolated elevation for the provided azimuth
+*/
+double DEMAlgorithms::getHorizon(const std::vector< std::pair<double,double> > &horizon, const double& azimuth)
+{
+	if (horizon.empty())
+		throw InvalidArgumentException("attempting to interpolate a horizon elevation from an empty set of horizon points", AT);
+	
+	//special (sick) case: only one value given for the horizon -> we consider this is a constant for all azimuths
+	if (horizon.size()==1)
+		return horizon.front().second;
+	
+	const std::vector< std::pair<double, double> >::const_iterator next = std::upper_bound(horizon.begin(), horizon.end(), make_pair(azimuth, 0.), sort_horizons()); //first element that is > azimuth
+	
+	double x1, y1, x2, y2;
+	if (next!=horizon.begin() && next!=horizon.end()) { //normal case
+		const size_t ii = next - horizon.begin();
+		x1 = horizon[ii-1].first;
+		y1 = horizon[ii-1].second;
+		x2 = horizon[ii].first;
+		y2 = horizon[ii].second;
+	} else {
+		x1 = horizon.back().first - 360.;
+		y1 = horizon.back().second;
+		x2 = horizon.front().first;
+		y2 = horizon.front().second;
 	}
+	
+	const double a = (y2 - y1) / (x2 - x1);
+	const double b = y2 - a * x2;
+	
+	return a*azimuth + b;
 }
 
 /**
diff --git a/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.h b/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.h
index e1594d74..417501f5 100644
--- a/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.h
+++ b/Source/meteoio/meteoio/dataClasses/DEMAlgorithms.h
@@ -33,8 +33,11 @@ class DEMAlgorithms {
 	public:
 		static Grid2DObject getHillshade(const DEMObject& dem, const double& elev, const double& azimuth);
 		static double getHorizon(const DEMObject& dem, const size_t& ix1, const size_t& iy1, const double& bearing);
-		static double getHorizon(const DEMObject& dem, const Coords& point, const double& bearing);
-		static void getHorizon(const DEMObject& dem, const Coords& point, const double& increment, std::vector< std::pair<double,double> >& horizon);
+		static double getHorizon(const DEMObject& dem, Coords point, const double& bearing);
+		static std::vector< std::pair<double,double> > getHorizonScan(const DEMObject& dem, Coords point, const double& increment);
+		static std::map< std::string, std::vector< std::pair<double,double> > > readHorizonScan(const std::string& where, const std::string& filename);
+		static double getHorizon(const std::vector< std::pair<double,double> > &horizon, const double& azimuth);
+		static void writeHorizons(const std::map< std::string, std::vector< std::pair<double,double> > >& horizon, const std::string& filename);
         static double getCellSkyViewFactor(const DEMObject& dem, const size_t& ii, const size_t& jj);
 
 	private:
diff --git a/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.cc b/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.cc
index 9189f4f3..0780a560 100644
--- a/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.cc
+++ b/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.cc
@@ -26,6 +26,13 @@
 
 namespace mio {
 
+AllSkyLWGenerator::AllSkyLWGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ, const Config &i_cfg)
+                   : TauCLDGenerator(vecArgs, i_algo, i_section, TZ, i_cfg), sun(), model(OMSTEDT)
+{ 
+	//TauCLDGenerator will do its own arguments parsing, then AllSkyLWGenerator
+	parse_args(vecArgs); 
+}
+
 void AllSkyLWGenerator::parse_args(const std::vector< std::pair<std::string, std::string> >& vecArgs)
 {
 	const std::string where( section+"::"+algo );
@@ -42,30 +49,18 @@ void AllSkyLWGenerator::parse_args(const std::vector< std::pair<std::string, std
 			else if (user_algo=="UNSWORTH") model = UNSWORTH;
 			else if (user_algo=="CRAWFORD") {
 				model = CRAWFORD;
-				cloudiness_model = TauCLDGenerator::CLF_CRAWFORD;
+				if (cloudiness_model==TauCLDGenerator::DEFAULT)
+					cloudiness_model = TauCLDGenerator::CLF_CRAWFORD;
 			} else if (user_algo=="LHOMME") {
 				model = LHOMME;
-				cloudiness_model = TauCLDGenerator::CLF_LHOMME;
+				if (cloudiness_model==TauCLDGenerator::DEFAULT)
+					cloudiness_model = TauCLDGenerator::CLF_LHOMME;
 			} else
 				throw InvalidArgumentException("Unknown parametrization \""+user_algo+"\" supplied for "+where, AT);
 
 			has_type = true;
-		} else if (vecArgs[ii].first=="CLOUDINESS_TYPE") {
-			user_cloudiness = IOUtils::strToUpper(vecArgs[ii].second);
-		} else if (vecArgs[ii].first=="USE_RSWR") {
-			IOUtils::parseArg(vecArgs[ii], where, use_rswr);
-		} else if (vecArgs[ii].first=="USE_RAD_THRESHOLD") {
-			IOUtils::parseArg(vecArgs[ii], where, use_rad_threshold);
 		}
 	}
-
-	if (!user_cloudiness.empty()) {
-		if (user_cloudiness=="LHOMME") cloudiness_model = TauCLDGenerator::CLF_LHOMME;
-		else if (user_cloudiness=="KASTEN") cloudiness_model = TauCLDGenerator::KASTEN;
-		else if (user_cloudiness=="CRAWFORD") cloudiness_model = TauCLDGenerator::CLF_CRAWFORD;
-		else
-			throw InvalidArgumentException("Unknown parametrization \""+user_cloudiness+"\" supplied for "+where, AT);
-	}
 	
 	if (!has_type) throw InvalidArgumentException("Please provide a TYPE for "+where, AT);
 }
@@ -98,7 +93,7 @@ bool AllSkyLWGenerator::generate(const size_t& param, MeteoData& md)
 			sun.setDate(julian_gmt, 0.);
 
 			bool is_night;
-			cloudiness = TauCLDGenerator::getCloudiness(cloudiness_model, md, use_rswr, use_rad_threshold, sun, is_night);
+			cloudiness = TauCLDGenerator::getCloudiness(md, sun, is_night);
 			if (cloudiness==IOUtils::nodata && !is_night) return false;
 
 			if (is_night) { //interpolate the cloudiness over the night
diff --git a/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.h b/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.h
index d4fd6bba..22c2cb9d 100644
--- a/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.h
+++ b/Source/meteoio/meteoio/dataGenerators/AllSkyLWGenerator.h
@@ -48,26 +48,14 @@ namespace mio {
  * for the Greenland Ice Sheet."</i> Global and Planetary change <b>9.1</b> (1994): 143-164.
  *      - UNSWORTH -- from Unsworth and Monteith, <i>"Long-wave radiation at the ground"</i>,
  * Q. J. R. Meteorolo. Soc., Vol. 101, 1975, pp 13-24 coupled with a clear sky emissivity following (Dilley, 1998).
- *  - CLOUDINESS_TYPE: normally, the cloudiness parametrization that might be needed to convert a clearness index (comparing the
- * measured ISWR to the potential ISWR, see TauCLDGenerator for more as well as the supported parametrizations) is given for 
- * each long wave parametrization but it is possible here to force it to a specific parametrization (default: the 
- * cloudiness parametrization that belongs to the long wave parametrization);
- *  - USE_RSWR: if set to TRUE, when no ISWR is available but RSWR and HS are available, a ground albedo is estimated
- * (either soil or snow albedo) and ISWR is then computed from RSWR. Unfortunatelly, this is not very precise... (thus default is false)
- *  - USE_RAD_THRESHOLD: when relying on measured ISWR to parametrize the cloudiness, there is a risk that the measuring station would
- * stand in a place where it is shaded by the surrounding terrain at some point during the day. This would lead to an overestimation 
- * of the cloudiness that is undesirable. In this case, it is possible to set USE_RAD_THRESHOLD to TRUE in order to interpolate the cloudiness
- * over all periods of low radiation measured ISWR. This is less performant that only considering the solar elevation but improves things
- * in this specific scenario.
+ *  - it also takes all the arguments of TauCLDGenerator, including the option to provide an Horizon file (see the ProcShade filter for the format).
  *
  * If no cloud transmissivity is provided in the data, it is calculated from the solar index (ratio of measured iswr to potential iswr, therefore using
- * the current location (lat, lon, altitude) and ISWR to parametrize the cloud cover). This relies on (Kasten and Czeplak, 1980)
- * except for Crawford and Lhomme that provide their own parametrizations.
- * The last evaluation of cloud transmissivity is used all along during the times when no ISWR is available if such ratio
- * is not too old (ie. no more than 1 day old).
- * If only RSWR is measured, the measured snow height is used to determine if there is snow on the ground or not.
- * In case of snow, a snow albedo of 0.85 is used while in the abscence of snow, a grass albedo of 0.23 is used
- * in order to compute ISWR from RSWR (please be aware that this very significantly degrades the performance of the parametrization).
+ * the current location (lat, lon, altitude) and ISWR to parametrize the cloud cover). This relies on (Kasten and Czeplak, 1980) by default 
+ * except for Crawford and Lhomme that provide their own parametrizations (it can be forced through the TauCLDGenerator options).
+ * The last evaluation of cloud transmissivity is used all along during the times when no ISWR is available (as it is night, or the station stands in 
+ * the shade or the ISWR measurement is missing) and the last valid value is not too old (ie. no more than 1 day old).
+ * 
  * Finally, it is recommended to also use a clear sky generator (declared after this one)
  * for the case of no available short wave measurement (by declaring the ClearSky generator \em after AllSky).
  * @code
@@ -84,11 +72,9 @@ namespace mio {
  * \image latex all_sky_ilwr_cmp.eps "Comparison between measured and parametrized ILWR at the Weissfluhjoch *WFJ station (2691m, Davos, Switzerland) for the 2010-08-01 – 2019-08-01 period" width=0.9\textwidth
  *
  */
-class AllSkyLWGenerator : public GeneratorAlgorithm {
+class AllSkyLWGenerator : public TauCLDGenerator {
 	public:
-		AllSkyLWGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ)
-		               : GeneratorAlgorithm(vecArgs, i_algo, i_section, TZ), sun(),
-		                 last_cloudiness(), model(OMSTEDT), cloudiness_model(TauCLDGenerator::KASTEN), use_rswr(false), use_rad_threshold(false) { parse_args(vecArgs); }
+		AllSkyLWGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ, const Config &i_cfg);
 		bool generate(const size_t& param, MeteoData& md);
 		bool create(const size_t& param, const size_t& ii_min, const size_t& ii_max, std::vector<MeteoData>& vecMeteo);
 	private:
@@ -104,10 +90,7 @@ class AllSkyLWGenerator : public GeneratorAlgorithm {
 		} parametrization;
 
 		SunObject sun;
-		std::map< std::string, std::pair<double, double> > last_cloudiness; //as < station_hash, <julian_gmt, cloudiness> >
 		parametrization model;
-		TauCLDGenerator::clf_parametrization cloudiness_model;
-		bool use_rswr, use_rad_threshold;
 };
 
 } //end namespace mio
diff --git a/Source/meteoio/meteoio/dataGenerators/GeneratorAlgorithms.cc b/Source/meteoio/meteoio/dataGenerators/GeneratorAlgorithms.cc
index c382a988..71d016a0 100644
--- a/Source/meteoio/meteoio/dataGenerators/GeneratorAlgorithms.cc
+++ b/Source/meteoio/meteoio/dataGenerators/GeneratorAlgorithms.cc
@@ -155,7 +155,7 @@ GeneratorAlgorithm* GeneratorAlgorithmFactory::getAlgorithm(const Config& cfg, c
 	} else if (algoname == "HUMIDITY"){
 		return new HumidityGenerator(vecArgs, i_algoname, i_section, TZ);
 	} else if (algoname == "TAU_CLD"){
-		return new TauCLDGenerator(vecArgs, i_algoname, i_section, TZ);
+		return new TauCLDGenerator(vecArgs, i_algoname, i_section, TZ, cfg);
 	} else if (algoname == "TS_OLWR"){
 		return new TsGenerator(vecArgs, i_algoname, i_section, TZ);
 	} else if (algoname == "ISWR_ALBEDO"){
@@ -163,7 +163,7 @@ GeneratorAlgorithm* GeneratorAlgorithmFactory::getAlgorithm(const Config& cfg, c
 	} else if (algoname == "CLEARSKY_LW"){
 		return new ClearSkyLWGenerator(vecArgs, i_algoname, i_section, TZ);
 	} else if (algoname == "ALLSKY_LW"){
-		return new AllSkyLWGenerator(vecArgs, i_algoname, i_section, TZ);
+		return new AllSkyLWGenerator(vecArgs, i_algoname, i_section, TZ, cfg);
 	} else if (algoname == "ALLSKY_SW"){
 		return new AllSkySWGenerator(vecArgs, i_algoname, i_section, TZ);
 	} else if (algoname == "CLEARSKY_SW"){
diff --git a/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.cc b/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.cc
index a7e55f93..a7acaf28 100644
--- a/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.cc
+++ b/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.cc
@@ -19,16 +19,21 @@
 
 #include <meteoio/dataGenerators/TauCLDGenerator.h>
 #include <meteoio/meteoLaws/Atmosphere.h>
+#include <meteoio/IOHandler.h>
+#include <meteoio/FileUtils.h>
+#include <meteoio/dataClasses/DEMAlgorithms.h>
 #include <algorithm>
 
 namespace mio {
 
-TauCLDGenerator::TauCLDGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ)
-                              : GeneratorAlgorithm(vecArgs, i_algo, i_section, TZ), last_cloudiness(), cloudiness_model(KASTEN), use_rswr(false), use_rad_threshold(false)
+TauCLDGenerator::TauCLDGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ, const Config &i_cfg)
+                              : GeneratorAlgorithm(vecArgs, i_algo, i_section, TZ), last_cloudiness(), masks(), horizons_outfile(), cfg(i_cfg), dem(), cloudiness_model(DEFAULT), use_rswr(false), use_rad_threshold(false), write_mask_out(), has_dem(false)
 {
 	const std::string where( section+"::"+algo );
+	bool from_dem=false, has_infile=false;
+	
 	for (size_t ii=0; ii<vecArgs.size(); ii++) {
-		if (vecArgs[ii].first=="TYPE") {
+		if (vecArgs[ii].first=="CLD_TYPE") {
 			const std::string user_algo( IOUtils::strToUpper(vecArgs[ii].second) );
 			
 			if (user_algo=="LHOMME") cloudiness_model = CLF_LHOMME;
@@ -43,9 +48,42 @@ TauCLDGenerator::TauCLDGenerator(const std::vector< std::pair<std::string, std::
 		if (vecArgs[ii].first=="USE_RAD_THRESHOLD") {
 			IOUtils::parseArg(vecArgs[ii], where, use_rad_threshold);
 		}
+		if (vecArgs[ii].first=="SHADE_FROM_DEM") {
+			IOUtils::parseArg(vecArgs[ii], where, from_dem);
+		}
+		if (vecArgs[ii].first=="INFILE") {
+			const std::string root_path( cfg.getConfigRootDir() );
+			//if this is a relative path, prefix the path with the current path
+			const std::string in_filename( vecArgs[ii].second );
+			const std::string prefix = ( FileUtils::isAbsolutePath(in_filename) )? "" : root_path+"/";
+			const std::string path( FileUtils::getPath(prefix+in_filename, true) );  //clean & resolve path
+			const std::string filename( path + "/" + FileUtils::getFilename(in_filename) );
+			masks = DEMAlgorithms::readHorizonScan(where, filename); //this mask is valid for ALL stations
+			has_infile = true;
+		} else if (vecArgs[ii].first=="OUTFILE") {
+			IOUtils::parseArg(vecArgs[ii], where, horizons_outfile);
+			write_mask_out = true;
+		}
+	}
+	
+	const bool use_default_horizon = !from_dem && !has_infile;
+	if (write_mask_out && use_default_horizon)
+		throw InvalidArgumentException("In "+where+", please provide an INFILE and/or a DEM to compute the shading and write it out!", AT);
+	
+	if (from_dem) {
+		IOHandler io(cfg);
+		dem.setUpdatePpt( DEMObject::NO_UPDATE ); //we only need the elevations
+		io.readDEM(dem);
+		has_dem = true;
 	}
 }
 
+TauCLDGenerator::~TauCLDGenerator()
+{
+	if (has_dem && !masks.empty() && !horizons_outfile.empty())
+		DEMAlgorithms::writeHorizons(masks, horizons_outfile);
+}
+
 /**
  * @brief Compute the clearness index from an atmospheric cloudiness value
  * @details
@@ -55,45 +93,75 @@ TauCLDGenerator::TauCLDGenerator(const std::vector< std::pair<std::string, std::
  * the same way as more traditional measurements (ie only ISWR provided) where it will be re-converted
  * to a cloudiness (thus falling abck to the same cloudiness as originally provided).
 
- * @param[in] clf_model cloudiness parametrization
  * @param[in] cloudiness cloudiness (between 0 and 1)
  * @return clearness index (between 0 and 1)
  */
-double TauCLDGenerator::getClearness(const clf_parametrization& clf_model, const double& cloudiness)
+double TauCLDGenerator::getClearness(const double& cloudiness) const
 {
-	if (clf_model==CLF_LHOMME) {
+	if (cloudiness_model==CLF_LHOMME) {
 		return Atmosphere::Lhomme_cloudiness( cloudiness/8. );
-	} else if (clf_model==KASTEN) {
+	} else if (cloudiness_model==KASTEN || cloudiness_model==DEFAULT) {
 		return Atmosphere::Kasten_cloudiness( cloudiness/8. );
-	} else if (clf_model==CLF_CRAWFORD) {
+	} else if (cloudiness_model==CLF_CRAWFORD) {
 		return Atmosphere::Lhomme_cloudiness( cloudiness/8. );
 	} else
 		return IOUtils::nodata; //this should never happen
 }
 
+std::vector< std::pair<double,double> > TauCLDGenerator::computeMask(const DEMObject& i_dem, const StationData& sd)
+{
+	//compute horizon by "angularResolution" increments
+	const double cellsize = i_dem.cellsize;
+	const double angularResolution = (cellsize<=20.)? 5. : (cellsize<=100.)? 10. : 20.; 
+	std::vector< std::pair<double,double> > o_mask( DEMAlgorithms::getHorizonScan(i_dem, sd.position, angularResolution) );
+	if (o_mask.empty()) throw InvalidArgumentException( "In Generator, could not compute mask from DEM '"+i_dem.llcorner.toString(Coords::LATLON)+"'", AT);
+
+	return o_mask;
+}
+
+
 /**
  * @brief Compute the atmospheric cloudiness from the available measurements
  * @details
  * The clearness index (ie the ratio of the incoming short wave radiation over the ground potential radiation, projected on the horizontal) 
  * is computed and used to evaluate the cloudiness, based on the chosen parametrization.
- * @param[in] clf_model cloudiness parametrization
  * @param[in] md MeteoData
- * @param[in] i_use_rswr if set to true, in case of no iswr measurements, a ground albedo is assumed and used to compute iswr. Based on HS, this albedo can either
  * be a soil ro a snow albedo
- * @param[in] i_use_rad_threshold use a radiation threshold to force resampling of cloudiness over prediods of low ISWR?
  * @param sun For better efficiency, the SunObject for this location (so it can be cached)
  * @param[out] is_night set to TRUE if it is night time
  * @return cloudiness (between 0 and 1)
  */
-double TauCLDGenerator::getCloudiness(const clf_parametrization& clf_model, const MeteoData& md, const bool& i_use_rswr, const bool& i_use_rad_threshold, SunObject& sun, bool &is_night)
+double TauCLDGenerator::getCloudiness(const MeteoData& md, SunObject& sun, bool &is_night)
 {
 	//we know that TA and RH are available, otherwise we would not get called
 	const double TA=md(MeteoData::TA), RH=md(MeteoData::RH), HS=md(MeteoData::HS), RSWR=md(MeteoData::RSWR);
 	double ISWR=md(MeteoData::ISWR);
-
+	
+	double sun_azi, sun_elev;
+	sun.position.getHorizontalCoordinates(sun_azi, sun_elev);
+	
+	//check if the station already has an associated mask, first by stationID then as wildcard
+	double mask_elev = 5.;
+	const std::string stationID( md.getStationID() );
+	std::map< std::string , std::vector< std::pair<double,double> > >::const_iterator mask = masks.find( stationID );
+	if (mask!=masks.end()) {
+		mask_elev = DEMAlgorithms::getHorizon(mask->second, sun_azi);
+	} else {
+		//now look for a wildcard fallback
+		mask = masks.find( "*" );
+		if (mask!=masks.end()) {
+			mask_elev = DEMAlgorithms::getHorizon(mask->second, sun_azi);
+		} else if (has_dem) {
+			masks[ stationID ] = computeMask(dem, md.meta);
+			mask = masks.find( stationID );
+			mask_elev = DEMAlgorithms::getHorizon(mask->second, sun_azi);
+		}
+	}
+	
 	//at sunrise or sunset, we might get very wrong results -> return nodata in order to use interpolation instead
 	//obviously, when it is really night neither can we compute anything here...
-	is_night = (sun.position.getSolarElevation() <= 5.);
+	//we use 5 degrees to represent very low elevation rays of light that would not work well in the horizontal sensors
+	is_night = (sun_elev <= mask_elev || sun_elev<5.);
 	if (is_night) return IOUtils::nodata;
 	
 	double albedo = .5;
@@ -104,7 +172,7 @@ double TauCLDGenerator::getCloudiness(const clf_parametrization& clf_model, cons
 			albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
 
 		if (ISWR==IOUtils::nodata) { //ISWR is missing, trying to compute it
-			if (!i_use_rswr) return IOUtils::nodata;
+			if (!use_rswr) return IOUtils::nodata;
 			if (RSWR!=IOUtils::nodata && HS!=IOUtils::nodata)
 				ISWR = RSWR / albedo;
 			else
@@ -117,7 +185,7 @@ double TauCLDGenerator::getCloudiness(const clf_parametrization& clf_model, cons
 	sun.getHorizontalRadiation(toa, direct, diffuse);
 	const double iswr_clear_sky = direct+diffuse;
 	
-	if (i_use_rad_threshold) {
+	if (use_rad_threshold) {
 		if (ISWR<Atmosphere::day_iswr_thresh || iswr_clear_sky<Atmosphere::day_iswr_thresh || iswr_clear_sky<ISWR) {
 			is_night = true;
 			return IOUtils::nodata;
@@ -125,15 +193,15 @@ double TauCLDGenerator::getCloudiness(const clf_parametrization& clf_model, cons
 	}
 
 	const double clearness_index = std::min(ISWR/iswr_clear_sky, 1.);
-	if (clf_model==CLF_LHOMME) {
+	if (cloudiness_model==CLF_LHOMME) {
 		const double clf = Atmosphere::Lhomme_cloudiness( clearness_index );
 		if (clf<0. || clf>1.) return IOUtils::nodata;
 		return clf;
-	} else if (clf_model==KASTEN) {
+	} else if (cloudiness_model==KASTEN || cloudiness_model==DEFAULT) {
 		const double clf = Atmosphere::Kasten_cloudiness( clearness_index );
 		if (clf<0. || clf>1.) return IOUtils::nodata;
 		return clf;
-	} else if (clf_model==CLF_CRAWFORD) {
+	} else if (cloudiness_model==CLF_CRAWFORD) {
 		const double clf = Atmosphere::Lhomme_cloudiness( clearness_index );
 		if (clf<0. || clf>1.) return IOUtils::nodata;
 		return clf;
@@ -149,7 +217,7 @@ bool TauCLDGenerator::generate(const size_t& param, MeteoData& md)
 		if (cld!=IOUtils::nodata) {
 			if (cld==9) cld=8.; //Synop sky obstructed from view -> fully cloudy
 			if (cld>8. || cld<0.) throw InvalidArgumentException("Cloud cover CLD should be between 0 and 8!", AT);
-			value = getClearness( cloudiness_model, cld );
+			value = getClearness( cld );
 			return true;
 		}
 
@@ -169,7 +237,7 @@ bool TauCLDGenerator::generate(const size_t& param, MeteoData& md)
 		sun.setDate(julian_gmt, 0.);
 
 		bool is_night;
-		double cloudiness = TauCLDGenerator::getCloudiness(cloudiness_model, md, use_rswr, use_rad_threshold, sun, is_night);
+		double cloudiness = TauCLDGenerator::getCloudiness(md, sun, is_night);
 		if (cloudiness==IOUtils::nodata && !is_night) return false;
 
 		if (is_night) { //interpolate the cloudiness over the night
diff --git a/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.h b/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.h
index eed981b6..ec74f63b 100644
--- a/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.h
+++ b/Source/meteoio/meteoio/dataGenerators/TauCLDGenerator.h
@@ -21,6 +21,7 @@
 
 #include <meteoio/dataGenerators/GeneratorAlgorithms.h>
 #include <meteoio/meteoLaws/Sun.h>
+#include <meteoio/dataClasses/DEMObject.h>
 
 #include <map>
 
@@ -41,8 +42,11 @@ namespace mio {
  * to parametrize the cloud cover). This relies on (Kasten and Czeplak, 1980).
  *
  * It takes the following (optional) argument:
- *    - TYPE: cloudiness model, either LHOMME, KASTEN or CRAWFORD (default: KASTEN, see AllSkyLWGenerator for the references of the papers. 
+ *    - CLD_TYPE: cloudiness model, either LHOMME, KASTEN or CRAWFORD (default: KASTEN, see AllSkyLWGenerator for the references of the papers. 
  * Please also note that CRAWFORD and LHOMME are exactly identical as the both simply consider that the cloudiness is <em>1-clearness_index</em>);
+ *    - SHADE_FROM_DEM: if set to true, the DEM defined in the [Input] section will be used to compute the shading;
+ *    - INFILE: a file containing the horizon for some or all stations (see the ProcShade filter for the format);
+ *    - OUTFILE: a file to write the computed horizons to. If some horizons have been read from INFILE, they will also be written out in OUTFILE;
  *    - USE_RSWR. If set to TRUE, when no ISWR is available but RSWR and HS are available, a ground albedo is estimated
  * (either soil or snow albedo) and ISWR is then computed from RSWR. Unfortunatelly, this is not very precise... (thus default is false)
  *    - USE_RAD_THRESHOLD: when relying on measured ISWR to parametrize the cloudiness, there is a risk that the measuring station would
@@ -50,6 +54,9 @@ namespace mio {
  * of the cloudiness that is undesirable. In this case, it is possible to set USE_RAD_THRESHOLD to TRUE in order to interpolate the cloudiness
  * over all periods of low radiation measured ISWR. This is less performant that only considering the solar elevation but improves things
  * in this specific scenario.
+ * 
+ * Please not that it is possible to combine SHADE_FROM_DEM and INFILE: in this case, stations that don't have any horizon in the provided
+ * INFILE will be computed from DEM. If none is available, a fixed 5 degrees threshold is used.
  *
  * @code
  * [Generators]
@@ -60,20 +67,30 @@ namespace mio {
 class TauCLDGenerator : public GeneratorAlgorithm {
 	public:
 		typedef enum CLF_PARAMETRIZATION {
+			DEFAULT,	//will be mapped to KASTEN
 			CLF_LHOMME,
 			KASTEN,
 			CLF_CRAWFORD
 		} clf_parametrization;
 
-		TauCLDGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ);
+		TauCLDGenerator(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& i_algo, const std::string& i_section, const double& TZ, const Config &i_cfg);
+		~TauCLDGenerator();
+		
 		bool generate(const size_t& param, MeteoData& md);
 		bool create(const size_t& param, const size_t& ii_min, const size_t& ii_max, std::vector<MeteoData>& vecMeteo);
-		static double getCloudiness(const clf_parametrization& clf_model, const MeteoData& md, const bool& i_use_rswr, const bool& i_use_rad_threshold, SunObject& sun, bool &is_night);
-		static double getClearness(const clf_parametrization& clf_model, const double& cloudiness);
-	private:
+	protected:
+		double getCloudiness(const MeteoData& md, SunObject& sun, bool &is_night);
+		double getClearness(const double& cloudiness) const;
+		static std::vector< std::pair<double,double> > computeMask(const DEMObject& i_dem, const StationData& sd);
+		
 		std::map< std::string, std::pair<double, double> > last_cloudiness; //as < station_hash, <julian_gmt, cloudiness> >
+		std::map< std::string , std::vector< std::pair<double,double> > > masks;
+		std::string horizons_outfile;
+		const Config &cfg;
+		DEMObject dem;
 		clf_parametrization cloudiness_model;
 		bool use_rswr, use_rad_threshold;
+		bool write_mask_out, has_dem;
 };
 
 } //end namespace mio
diff --git a/Source/meteoio/meteoio/meteoFilters/ProcShade.cc b/Source/meteoio/meteoio/meteoFilters/ProcShade.cc
index 136ebe40..f2c2e467 100644
--- a/Source/meteoio/meteoio/meteoFilters/ProcShade.cc
+++ b/Source/meteoio/meteoio/meteoFilters/ProcShade.cc
@@ -16,54 +16,50 @@
     You should have received a copy of the GNU Lesser General Public License
     along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
 */
-#include <fstream>
-#include <sstream>
-#include <cerrno>
-#include <cstring>
-#include <algorithm>
 
 #include <meteoio/meteoFilters/ProcShade.h>
 #include <meteoio/meteoLaws/Sun.h>
 #include <meteoio/IOHandler.h>
+#include <meteoio/IOUtils.h>
 #include <meteoio/FileUtils.h>
 #include <meteoio/dataClasses/DEMAlgorithms.h>
 
 using namespace std;
 
 namespace mio {
-//custom function for sorting cache_meteo_files
-struct sort_pred {
-	bool operator()(const std::pair<double,double> &left, const std::pair<double,double> &right) {
-		if (left.first < right.first) return true; else return false;
-	}
-};
 
 const double ProcShade::diffuse_thresh = 15.; //below this threshold, not correction is performed since it will only be diffuse
 
 ProcShade::ProcShade(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& name, const Config& i_cfg)
-        : ProcessingBlock(vecArgs, name, i_cfg), cfg(i_cfg), dem(), masks(), write_mask_out(false)
+        : ProcessingBlock(vecArgs, name, i_cfg), cfg(i_cfg), dem(), masks(), horizons_outfile(), has_dem(false), write_mask_out(false)
 {
 	parse_args(vecArgs);
 	properties.stage = ProcessingProperties::first; //for the rest: default values
 }
 
+ProcShade::~ProcShade()
+{
+	if (has_dem && !masks.empty() && !horizons_outfile.empty())
+		DEMAlgorithms::writeHorizons(masks, horizons_outfile);
+}
+
 void ProcShade::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
                         std::vector<MeteoData>& ovec)
 {
 	ovec = ivec;
 	if (ovec.empty()) return;
 	
-	const std::string stationHash( ovec[0].meta.getHash() );
 	SunObject Sun;
 	
-	//check if the station already has an associated mask, first as wildcard then by station hash
-	std::map< std::string , std::vector< std::pair<double,double> > >::iterator mask = masks.find( "*" );
+	//check if the station already has an associated mask, first by stationID then as wildcard
+	const std::string stationID( ovec[0].getStationID() );
+	std::map< std::string , std::vector< std::pair<double,double> > >::const_iterator mask = masks.find( stationID );
 	if (mask==masks.end()) {
-		//now look for our specific station hash
-		mask = masks.find( stationHash );
+		//now look for a wildcard fallback
+		mask = masks.find( "*" );
 		if (mask==masks.end()) {
-			masks[ stationHash ] = computeMask(dem, ovec[0].meta, write_mask_out);
-			mask = masks.find( stationHash);
+			masks[ stationID ] = computeMask(dem, ovec[0].meta);
+			mask = masks.find( stationID);
 		}
 	}
 	
@@ -80,7 +76,7 @@ void ProcShade::process(const unsigned int& param, const std::vector<MeteoData>&
 		double sun_azi, sun_elev;
 		Sun.position.getHorizontalCoordinates(sun_azi, sun_elev);
 		
-		const double mask_elev = getMaskElevation(mask->second, sun_azi);
+		const double mask_elev = DEMAlgorithms::getHorizon(mask->second, sun_azi);
 		if (mask_elev>0 && mask_elev>sun_elev) { //the point is in the shade
 			const double TA=ovec[ii](MeteoData::TA), RH=ovec[ii](MeteoData::RH), HS=ovec[ii](MeteoData::HS), RSWR=ovec[ii](MeteoData::RSWR);
 			double ISWR=ovec[ii](MeteoData::ISWR);
@@ -114,104 +110,14 @@ void ProcShade::process(const unsigned int& param, const std::vector<MeteoData>&
 	}
 }
 
-//linear interpolation between the available points
-double ProcShade::getMaskElevation(const std::vector< std::pair<double,double> > &mask, const double& azimuth) const
-{
-	const std::vector< std::pair<double, double> >::const_iterator next = std::upper_bound(mask.begin(), mask.end(), make_pair(azimuth, 0.), sort_pred()); //first element that is > azimuth
-	
-	double x1, y1, x2, y2;
-	if (next!=mask.begin() && next!=mask.end()) { //normal case
-		const size_t ii = next - mask.begin();
-		x1 = mask[ii-1].first;
-		y1 = mask[ii-1].second;
-		x2 = mask[ii].first;
-		y2 = mask[ii].second;
-	} else {
-		x1 = mask.back().first - 360.;
-		y1 = mask.back().second;
-		x2 = mask.front().first;
-		y2 = mask.front().second;
-	}
-	
-	const double a = (y2 - y1) / (x2 - x1);
-	const double b = y2 - a * x2;
-	
-	return a*azimuth + b;
-}
-
-std::vector< std::pair<double,double> > ProcShade::readMask(const std::string& filter, const std::string& filename)
+std::vector< std::pair<double,double> > ProcShade::computeMask(const DEMObject& i_dem, const StationData& sd)
 {
-	std::vector< std::pair<double,double> > o_mask;
-	std::ifstream fin(filename.c_str());
-	if (fin.fail()) {
-		std::ostringstream ss;
-		ss << "Filter " << filter << ": " << "error opening file \"" << filename << "\", possible reason: " << std::strerror(errno);
-		throw AccessException(ss.str(), AT);
-	}
-
-	const char eoln = FileUtils::getEoln(fin); //get the end of line character for the file
-
-	try {
-		size_t lcount=0;
-		double azimuth, value;
-		std::string line;
-		do {
-			lcount++;
-			getline(fin, line, eoln); //read complete line
-			IOUtils::stripComments(line);
-			IOUtils::trim(line);
-			if (line.empty()) continue;
-
-			std::istringstream iss(line);
-			iss.setf(std::ios::fixed);
-			iss.precision(std::numeric_limits<double>::digits10);
-			iss >> std::skipws >> azimuth;
-			if ( !iss || azimuth<0. || azimuth>360.) {
-				std::ostringstream ss;
-				ss << "Invalid azimuth in file " << filename << " at line " << lcount;
-				throw InvalidArgumentException(ss.str(), AT);
-			}
-			iss >> std::skipws >> value;
-			if ( !iss ){
-				std::ostringstream ss;
-				ss << "Invalid value in file " << filename << " at line " << lcount;
-				throw InvalidArgumentException(ss.str(), AT);
-			}
-
-			o_mask.push_back( make_pair(azimuth, value) );
-		} while (!fin.eof());
-		fin.close();
-	} catch (const std::exception&){
-		if (fin.is_open()) {//close fin if open
-			fin.close();
-		}
-		throw;
-	}
-	
-	if (o_mask.empty()) throw InvalidArgumentException("In filter 'SHADE', no valid mask found in file '"+filename+"'", AT);
-	std::sort(o_mask.begin(), o_mask.end(), sort_pred());
-	return o_mask;
-}
-
-std::vector< std::pair<double,double> > ProcShade::computeMask(const DEMObject& i_dem, const StationData& sd, const bool& dump_mask)
-{
-	std::vector< std::pair<double,double> > o_mask;
-	Coords position( sd.position );
-	if (!i_dem.gridify(position)) {
-		const string msg = "In filter 'SHADE', station '"+sd.stationID+"' "+position.toString(Coords::LATLON)+" is not included in the DEM "+i_dem.llcorner.toString(Coords::LATLON);
-		throw NoDataException(msg, AT);
-	}
-	
-	DEMAlgorithms::getHorizon(i_dem, position, 10., o_mask); //by 10deg increments
+	//compute horizon by "angularResolution" increments
+	const double cellsize = i_dem.cellsize;
+	const double angularResolution = (cellsize<=20.)? 5. : (cellsize<=100.)? 10. : 20.;
+	std::vector< std::pair<double,double> > o_mask( DEMAlgorithms::getHorizonScan(i_dem, sd.position, angularResolution) );
 	if (o_mask.empty()) throw InvalidArgumentException( "In filter 'SHADE', could not compute mask from DEM '"+i_dem.llcorner.toString(Coords::LATLON)+"'", AT);
 
-	if (dump_mask) {
-		std::cout << "Horizon mask for station '" << sd.stationID << "'\n";
-		for (size_t ii=0; ii<o_mask.size(); ii++)
-			std::cout << o_mask[ii].first << " " << o_mask[ii].second << "\n";
-		std::cout << "\n";
-	}
-
 	return o_mask;
 }
 
@@ -221,17 +127,17 @@ void ProcShade::parse_args(const std::vector< std::pair<std::string, std::string
 	bool from_dem=true;
 
 	for (size_t ii=0; ii<vecArgs.size(); ii++) {
-		if (vecArgs[ii].first=="FILE") {
+		if (vecArgs[ii].first=="INFILE") {
 			const std::string root_path( cfg.getConfigRootDir() );
 			//if this is a relative path, prefix the path with the current path
 			const std::string in_filename( vecArgs[ii].second );
 			const std::string prefix = ( FileUtils::isAbsolutePath(in_filename) )? "" : root_path+"/";
 			const std::string path( FileUtils::getPath(prefix+in_filename, true) );  //clean & resolve path
 			const std::string filename( path + "/" + FileUtils::getFilename(in_filename) );
-			masks["*"] = readMask(getName(), filename); //this mask is valid for ALL stations
+			masks = DEMAlgorithms::readHorizonScan(getName(), filename); //this mask is valid for ALL stations
 			from_dem = false;
-		} else if (vecArgs[ii].first=="DUMP_MASK") {
-			IOUtils::parseArg(vecArgs[ii], where, write_mask_out);
+		} else if (vecArgs[ii].first=="OUTFILE") {
+			IOUtils::parseArg(vecArgs[ii], where, horizons_outfile);
 		}
 	}
 
@@ -239,6 +145,7 @@ void ProcShade::parse_args(const std::vector< std::pair<std::string, std::string
 		IOHandler io(cfg);
 		dem.setUpdatePpt( DEMObject::NO_UPDATE ); //we only need the elevations
 		io.readDEM(dem);
+		has_dem=true;
 	}
 }
 
diff --git a/Source/meteoio/meteoio/meteoFilters/ProcShade.h b/Source/meteoio/meteoio/meteoFilters/ProcShade.h
index d941a96c..27152f2d 100644
--- a/Source/meteoio/meteoio/meteoFilters/ProcShade.h
+++ b/Source/meteoio/meteoio/meteoFilters/ProcShade.h
@@ -28,7 +28,7 @@
 namespace mio {
 
 /**
- * @class  ProcShade
+ * @class ProcShade
  * @ingroup processing
  * @brief Apply a shading mask to the Incoming or Reflected Short Wave Radiation
  * @details
@@ -37,23 +37,28 @@ namespace mio {
  * This mask will be linearly interpolated between the provided points in order to be applied to the true sun position.
  * 
  * When providing the shading mask in a separate file, the same mask will be applied to all stations. It simply need to
- * contain the horizon elevation (in deg.) as a function of azimuth (in deg.):
+ * contain the horizon elevation (in deg.) as a function of azimuth (in deg.), prefixed by either a stationID or "*" as fallback 
+ * wildcard for all stationIDs:
  * @code
- * 0	5
- * 15	25
- * 45	12
- * 180	30
- * 270	20
+ * DAV1	0	5
+ * DAV1	15	25
+ * DAV1	45	12
+ * DAV1	180	30
+ * DAV1	270	20
+ * SLF2	0	7.5
+ * SLF2	180	10
+ * *	0	5
+ * *	180	5
  * @endcode
  *
  * Therefore, the filter supports the following arguments:
- *  - FILE: file (and path) to read the mask from (optional);
- *  - DUMP_MASK: if set to TRUE, each mask will be printed on the screen as it is computed (optional).
+ *  - INFILE: file (and path) to read the mask from (optional);
+ *  - OUTFILE: file (and path) to write the mask to (optional, this is usefull if some masks have been computed from the DEM).
  * 
  * Then the filter is declared with the file name containing the horizon mask as argument:
  * @code
  * ISWR::filter1    = SHADE
- * ISWR::arg1::file = ../input/iswr_mask.dat
+ * ISWR::arg1::infile = ../input/iswr_mask.dat
  * @endcode
  *
  * If no arguments are provided, then it will compute the mask from the Digital Elevation Model. In such as case, 
@@ -64,21 +69,20 @@ namespace mio {
 class ProcShade : public ProcessingBlock {
 	public:
 		ProcShade(const std::vector< std::pair<std::string, std::string> >& vecArgs, const std::string& name, const Config &i_cfg);
+		~ProcShade();
 
 		virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
 		                     std::vector<MeteoData>& ovec);
-		
-		static std::vector< std::pair<double,double> > computeMask(const DEMObject& i_dem, const StationData& sd, const bool& dump_mask=false);
 
 	private:
-		static std::vector< std::pair<double,double> > readMask(const std::string& filter, const std::string& filename);
 		void parse_args(const std::vector< std::pair<std::string, std::string> >& vecArgs);
-		double getMaskElevation(const std::vector< std::pair<double,double> > &mask, const double& azimuth) const;
+		static std::vector< std::pair<double,double> > computeMask(const DEMObject& i_dem, const StationData& sd);
 
 		const Config &cfg;
 		DEMObject dem;
 		std::map< std::string , std::vector< std::pair<double,double> > > masks;
-		bool write_mask_out;
+		std::string horizons_outfile;
+		bool has_dem, write_mask_out;
 
 		static const double diffuse_thresh;
 };
diff --git a/Source/meteoio/meteoio/plugins/CMakeLists.txt b/Source/meteoio/meteoio/plugins/CMakeLists.txt
index 54e114ea..bced24f8 100644
--- a/Source/meteoio/meteoio/plugins/CMakeLists.txt
+++ b/Source/meteoio/meteoio/plugins/CMakeLists.txt
@@ -68,6 +68,21 @@ IF(PLUGIN_IMISIO)
 	SET(plugins_sources ${plugins_sources} plugins/ImisIO.cc)
 ENDIF(PLUGIN_IMISIO)
 
+IF(PLUGIN_METEOBLUE)
+	FIND_PACKAGE(CURL REQUIRED)
+	INCLUDE_DIRECTORIES(SYSTEM ${CURL_INCLUDE_DIRS})
+	SET(plugin_libs ${plugin_libs} ${CURL_LIBRARIES})
+	SET(plugins_sources ${plugins_sources} plugins/MeteoBlue.cc)
+ENDIF(PLUGIN_METEOBLUE)
+
+IF(PLUGIN_NETCDFIO)
+	FIND_PACKAGE(NetCDF REQUIRED)
+	INCLUDE_DIRECTORIES(SYSTEM ${NETCDF_INCLUDE_DIR})
+	SET(plugin_libs ${plugin_libs} ${NETCDF_LIBRARIES})
+	SET(plugins_sources ${plugins_sources} plugins/libncpp.cc)
+	SET(plugins_sources ${plugins_sources} plugins/NetCDFIO.cc)
+ENDIF(PLUGIN_NETCDFIO)
+
 IF(PLUGIN_OSHDIO)
 	FIND_PACKAGE(Matio REQUIRED)
 	INCLUDE_DIRECTORIES(SYSTEM ${Matio_INCLUDE_DIR})
@@ -104,21 +119,6 @@ IF(PLUGIN_SASEIO)
 	SET(plugins_sources ${plugins_sources} plugins/SASEIO.cc)
 ENDIF(PLUGIN_SASEIO)
 
-IF(PLUGIN_METEOBLUE)
-	FIND_PACKAGE(CURL REQUIRED)
-	INCLUDE_DIRECTORIES(SYSTEM ${CURL_INCLUDE_DIRS})
-	SET(plugin_libs ${plugin_libs} ${CURL_LIBRARIES})
-	SET(plugins_sources ${plugins_sources} plugins/MeteoBlue.cc)
-ENDIF(PLUGIN_METEOBLUE)
-
-IF(PLUGIN_NETCDFIO)
-	FIND_PACKAGE(NetCDF REQUIRED)
-	INCLUDE_DIRECTORIES(SYSTEM ${NETCDF_INCLUDE_DIR})
-	SET(plugin_libs ${plugin_libs} ${NETCDF_LIBRARIES})
-	SET(plugins_sources ${plugins_sources} plugins/libncpp.cc)
-	SET(plugins_sources ${plugins_sources} plugins/NetCDFIO.cc)
-ENDIF(PLUGIN_NETCDFIO)
-
 IF(PLUGIN_SMETIO)
 	SET(plugins_sources ${plugins_sources} plugins/SMETIO.cc)
 ENDIF(PLUGIN_SMETIO)
@@ -127,6 +127,13 @@ IF(PLUGIN_SNIO)
 	SET(plugins_sources ${plugins_sources} plugins/SNIO.cc)
 ENDIF(PLUGIN_SNIO)
 
+IF(PLUGIN_WWCSIO)
+	FIND_PACKAGE(MySQL REQUIRED)
+	INCLUDE_DIRECTORIES(SYSTEM ${MYSQL_INCLUDE_DIR})
+	SET(plugin_libs ${plugin_libs} ${MYSQL_LIBRARY})
+	SET(plugins_sources ${plugins_sources} plugins/WWCSIO.cc)
+ENDIF(PLUGIN_WWCSIO)
+
 IF(PLUGIN_ZRXPIO)
 	SET(plugins_sources ${plugins_sources} plugins/ZRXPIO.cc)
 ENDIF(PLUGIN_ZRXPIO)
diff --git a/Source/meteoio/meteoio/plugins/WWCSIO.cc b/Source/meteoio/meteoio/plugins/WWCSIO.cc
new file mode 100644
index 00000000..a6c71887
--- /dev/null
+++ b/Source/meteoio/meteoio/plugins/WWCSIO.cc
@@ -0,0 +1,350 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+/***********************************************************************************/
+/*  Copyright 2014 Snow and Avalanche Study Establishment    SASE-CHANDIGARH       */
+/***********************************************************************************/
+/* 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/>.
+*/
+#include <meteoio/plugins/WWCSIO.h>
+
+#ifdef _WIN32
+	#include <winsock.h>
+#endif // _WIN32
+
+#include <mysql.h>
+#include <stdio.h>
+
+#include <algorithm>
+
+using namespace std;
+
+namespace mio {
+/**
+* @page wwcs WWCSIO
+* @section WWCS_format Format
+* This is the plugin required to get meteorological data from the WWCS MySQL database.
+*
+* @section WWCS_units Units
+* The units are assumed to be the following:
+* - __temperatures__ in celsius
+* - __relative humidity__ in %
+* - __wind speed__ in m/s
+* - __precipitations__ in mm/h
+* - __radiation__ in W/m²
+*
+* @section WWCS_keywords Keywords
+* This plugin uses the following keywords:
+* - COORDSYS: coordinate system (see Coords); [Input] and [Output] section
+* - COORDPARAM: extra coordinates parameters (see Coords); [Input] and [Output] section
+* - WWCS_HOST: MySQL Host Name (e.g. localhost or 191.168.145.20); [Input] section
+* - WWCS_DB: MySQL Database (e.g. snowpack); [Input] section
+* - WWCS_USER: MySQL User Name (e.g. root); [Input] section
+* - WWCS_PASS: MySQL password; [Input] section
+* - TIME_ZONE: For [Input] and [Output] sections
+* - STATION#: station code for the given number #; [Input] section
+*/
+
+const double WWCSIO::plugin_nodata = -999.; //plugin specific nodata value. It can also be read by the plugin (depending on what is appropriate)
+const string WWCSIO::MySQLQueryStationMetaData = "SELECT stationName, latitude, longitude, altitude, slope, azimuth FROM metadata WHERE StationID='STATION_ID' ";
+const string WWCSIO::MySQLQueryMeteoData = "SELECT timestamp, ta, rh, p, logger_ta, logger_rh FROM meteoseries WHERE stationID='STATION_NUMBER' AND TimeStamp>='START_DATE' AND TimeStamp<='END_DATE' ORDER BY timestamp ASC";
+
+WWCSIO::WWCSIO(const std::string& configfile)
+        : cfg(configfile), vecStationIDs(), vecStationMetaData(),
+          mysqlhost(), mysqldb(), mysqluser(), mysqlpass(),
+          coordin(), coordinparam(), coordout(), coordoutparam(),
+          in_dflt_TZ(5.), out_dflt_TZ(5.)
+{
+	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
+	readConfig();
+}
+
+WWCSIO::WWCSIO(const Config& cfgreader)
+        : cfg(cfgreader), vecStationIDs(), vecStationMetaData(),
+          mysqlhost(), mysqldb(), mysqluser(), mysqlpass(),
+          coordin(), coordinparam(), coordout(), coordoutparam(),
+          in_dflt_TZ(5.), out_dflt_TZ(5.)
+{
+	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
+	readConfig();
+}
+
+void WWCSIO::readConfig()
+{
+	cfg.getValue("WWCS_HOST", "Input", mysqlhost);
+	cfg.getValue("WWCS_DB", "Input", mysqldb);
+	cfg.getValue("WWCS_USER", "Input", mysqluser);
+	cfg.getValue("WWCS_PASS", "Input", mysqlpass);
+	cfg.getValue("TIME_ZONE","Input", in_dflt_TZ, IOUtils::nothrow);
+	cfg.getValue("TIME_ZONE","Output", out_dflt_TZ, IOUtils::nothrow);
+}
+
+void WWCSIO::readStationIDs(std::vector<std::string>& vecStationID) const
+{
+	vecStationID.clear();
+	cfg.getValues("STATION", "INPUT", vecStationID);
+
+	if (vecStationID.empty()) {
+		cerr << "\tNo stations specified for WWCSIO... is this what you want?\n";
+	}
+}
+
+void WWCSIO::readStationMetaData()
+{
+	vecStationMetaData.clear();
+	std::vector<std::string> vecStationID;
+	readStationIDs( vecStationID );
+
+	for (size_t ii=0; ii<vecStationID.size(); ii++) {
+		// Retrieve the station IDs - this only needs to be done once per instance
+		std::string stat_abk, stao_nr;
+		parseStationID(vecStationID[ii], stat_abk, stao_nr);
+
+		// Retrieve the station meta data - this only needs to be done once per instance
+		std::vector<std::string> stationMetaData;
+		getStationMetaData(stat_abk, stao_nr, MySQLQueryStationMetaData, stationMetaData);
+
+		//SELECT stationName, latitude, longitude, altitude, slope, azimuth FROM metadata WHERE StationID='STATION_ID'
+		std::string stao_name;
+		double lat, longi, alt, slope, azi;
+		if (!IOUtils::convertString(stao_name, stationMetaData.at(0)))
+			throw ConversionFailedException("Invalid station name for station "+vecStationID[ii]+": "+stationMetaData.at(0), AT);
+		if (!IOUtils::convertString(lat, stationMetaData.at(1), std::dec))
+			throw ConversionFailedException("Invalid latitude for station "+vecStationID[ii]+": "+stationMetaData.at(1), AT);
+		if (!IOUtils::convertString(longi, stationMetaData.at(2), std::dec))
+			throw ConversionFailedException("Invalid longitude for station "+vecStationID[ii]+": "+stationMetaData.at(2), AT);
+		if (!IOUtils::convertString(alt, stationMetaData.at(3), std::dec))
+			throw ConversionFailedException("Invalid altitude for station "+vecStationID[ii]+": "+stationMetaData.at(3), AT);
+		if (!IOUtils::convertString(slope, stationMetaData.at(4), std::dec))
+			throw ConversionFailedException("Invalid slope for station "+vecStationID[ii]+": "+stationMetaData.at(3), AT);
+		if (!IOUtils::convertString(azi, stationMetaData.at(5), std::dec))
+			throw ConversionFailedException("Invalid azimuth for station "+vecStationID[ii]+": "+stationMetaData.at(3), AT);
+
+		Coords location(coordin,coordinparam);
+		location.setLatLon(lat, longi, alt);
+		const std::string station_name = (!stao_name.empty())? vecStationID[ii] + ":" + stao_name : vecStationID[ii];
+		StationData sd(location, vecStationID[ii], station_name);
+		sd.setSlope(slope, azi);
+		vecStationMetaData.push_back( sd );
+	}
+
+}
+
+void WWCSIO::getStationMetaData(const std::string& stat_abk, const std::string& stao_nr,
+                                 const std::string& sqlQuery, std::vector<std::string>& vecMetaData)
+{
+	vecMetaData.clear();
+
+	MYSQL *conn = mysql_init(nullptr);
+	if (!mysql_real_connect(conn, mysqlhost.c_str(), mysqluser.c_str(), mysqlpass.c_str(), mysqldb.c_str(), 0, nullptr, 0)) {
+		throw IOException("Could not initiate connection to Mysql server "+mysqlhost, AT);
+	}
+
+	std::string Query1( sqlQuery );
+	const std::string str1( "STATION_ID" );
+	Query1.replace(Query1.find(str1), str1.length(), stat_abk);// set 1st variable's value
+	if (mysql_query(conn, Query1.c_str())) {
+		throw IOException("Query1 \""+Query1+"\" failed to execute", AT);
+	}
+
+	MYSQL_RES *res = mysql_use_result(conn);
+	const unsigned int column_no = mysql_num_fields(res);
+	std::string tmp_str;
+	MYSQL_ROW row;
+	while ( ( row = mysql_fetch_row(res) ) != nullptr ) {
+		for (unsigned int ii=0; ii<column_no; ii++) {
+			IOUtils::convertString(tmp_str, row[ii]);
+			vecMetaData.push_back(tmp_str);
+		}
+	}
+	mysql_free_result(res);
+	mysql_close(conn);
+
+	const size_t nr_metadata = vecMetaData.size();
+	if (nr_metadata==0)
+		throw NoDataException("Station " + stat_abk+stao_nr + " not found in the database", AT);
+	if (nr_metadata<4)
+		throw ConversionFailedException("Error while converting station meta data for station "+stat_abk+stao_nr, AT);
+}
+
+void WWCSIO::parseStationID(const std::string& stationID, std::string& stat_abk, std::string& stao_nr)
+{
+	stat_abk = stationID.substr(0, stationID.length()-1); //The station name: e.g. KLO
+	stao_nr = stationID.substr(stationID.length()-1, 1); //The station number: e.g. 2
+	if (!std::isdigit(stao_nr[0])) {
+		//the station is one of these non-imis stations that don't contain a number...
+		stat_abk = stationID;
+		stao_nr = "0";
+	}
+}
+
+void WWCSIO::readStationData(const Date&, std::vector<StationData>& vecStation)
+{
+	vecStation.clear();
+	readStationMetaData(); //reads all the station meta data into the vecStationMetaData (member vector)
+	vecStation = vecStationMetaData; //vecStationMetaData is a global vector holding all meta data
+}
+
+void WWCSIO::readMeteoData(const Date& dateStart , const Date& dateEnd,
+                            std::vector<std::vector<MeteoData> >& vecMeteo)
+{
+	if (vecStationMetaData.empty()) readStationMetaData();
+
+	vecMeteo.clear();
+	vecMeteo.insert(vecMeteo.begin(), vecStationMetaData.size(), vector<MeteoData>());
+
+	for (size_t ii=0; ii<vecStationMetaData.size(); ii++) { //loop through relevant stations
+		readData(dateStart, dateEnd, vecMeteo, ii, vecStationMetaData);
+	}
+}
+
+//read meteo data for one station
+void WWCSIO::readData(const Date& dateStart, const Date& dateEnd, std::vector< std::vector<MeteoData> >& vecMeteo,
+                       const size_t& stationindex, const std::vector<StationData>& vecMeta) const
+{
+	vecMeteo.at(stationindex).clear();
+
+	Date dateS(dateStart), dateE(dateEnd);
+	dateS.setTimeZone(in_dflt_TZ);
+	dateE.setTimeZone(in_dflt_TZ);
+
+	std::string stat_abk, stao_nr;
+	std::vector<std::string> vecHTS1;
+	std::vector< std::vector<std::string> > vecResult;
+	parseStationID(vecMeta.at(stationindex).getStationID(), stat_abk, stao_nr);
+
+	getStationData(stat_abk, stao_nr, dateS, dateE, vecHTS1, vecResult);
+	MeteoData tmpmd;
+	tmpmd.meta = vecMeta.at(stationindex);
+
+	for (size_t ii=0; ii<vecResult.size(); ii++){
+		parseDataSet(vecResult[ii], tmpmd);
+		convertUnits(tmpmd);
+		vecMeteo.at(stationindex).push_back( tmpmd ); //Now insert tmpmd
+	}
+}
+
+void WWCSIO::convertUnits(MeteoData& meteo)
+{
+	meteo.standardizeNodata(plugin_nodata);
+
+	//converts C to Kelvin, converts RH to [0,1]
+	double& ta = meteo(MeteoData::TA);
+	ta = IOUtils::C_TO_K(ta);
+
+	double& tss = meteo(MeteoData::TSS);
+	tss = IOUtils::C_TO_K(tss);
+
+	double& rh = meteo(MeteoData::RH);
+	if (rh != IOUtils::nodata)
+	rh /= 100.;
+
+	/*double& hs = meteo(MeteoData::HS);
+	if (hs != IOUtils::nodata)
+	hs /= 100.0;
+
+	double& rswr = meteo(MeteoData::RSWR);
+	if (rswr != IOUtils::nodata)
+	rswr /= 100.0;*/
+}
+
+void WWCSIO::parseDataSet(const std::vector<std::string>& i_meteo, MeteoData& md) const
+{
+	const std::string statID( md.meta.getStationID() );
+
+	if (!IOUtils::convertString(md.date, i_meteo.at(0), in_dflt_TZ))
+		throw ConversionFailedException("Invalid timestamp for station "+statID+": "+i_meteo.at(0), AT);
+	if (!IOUtils::convertString(md(MeteoData::TA), i_meteo.at(1)))
+		throw ConversionFailedException("Invalid TA for station "+statID+": "+i_meteo.at(1), AT);
+	if (!IOUtils::convertString(md(MeteoData::RH), i_meteo.at(2)))
+		throw ConversionFailedException("Invalid RH for station "+statID+": "+i_meteo.at(2), AT);
+	if (!IOUtils::convertString(md(MeteoData::VW), i_meteo.at(3)))
+		throw ConversionFailedException("Invalid VW for station "+statID+": "+i_meteo.at(3), AT);
+	if (!IOUtils::convertString(md(MeteoData::HS), i_meteo.at(4)))
+		throw ConversionFailedException("Invalid HS for station "+statID+": "+i_meteo.at(4), AT);
+	if (!IOUtils::convertString(md(MeteoData::TSS), i_meteo.at(5)))
+		throw ConversionFailedException("Invalid TSS for station "+statID+": "+i_meteo.at(5), AT);
+	if (!IOUtils::convertString(md(MeteoData::RSWR), i_meteo.at(6)))
+		throw ConversionFailedException("Invalid RSWR for station "+statID+": "+i_meteo.at(6), AT);
+	if (!IOUtils::convertString(md(MeteoData::ISWR), i_meteo.at(7)))
+		throw ConversionFailedException("Invalid ISWR for station "+statID+": "+i_meteo.at(7), AT);
+	if (!IOUtils::convertString(md(MeteoData::ILWR), i_meteo.at(8)))
+		throw ConversionFailedException("Invalid ILWR for station "+statID+": "+i_meteo.at(8), AT);
+	if (!IOUtils::convertString(md(MeteoData::PSUM), i_meteo.at(9)))
+		throw ConversionFailedException("Invalid PSUM for station "+statID+": "+i_meteo.at(9), AT);
+	if (!IOUtils::convertString(md(MeteoData::P), i_meteo.at(10)))
+		throw ConversionFailedException("Invalid P for station "+statID+": "+i_meteo.at(10), AT);
+}
+
+bool WWCSIO::getStationData(const std::string& stat_abk, const std::string& stao_nr,
+                             const Date& dateS, const Date& dateE,
+                             const std::vector<std::string>& vecHTS1,
+                             std::vector< std::vector<std::string> >& vecMeteoData) const
+{
+	vecMeteoData.clear();
+	bool fullStation = true;
+
+	// Creating MySQL TimeStamps
+	std::string sDate( dateS.toString(Date::ISO) );
+	std::replace( sDate.begin(), sDate.end(), 'T', ' ');
+	std::string eDate( dateE.toString(Date::ISO) );
+	std::replace( eDate.begin(), eDate.end(), 'T', ' ');
+
+	//SELECT timestamp, ta, rh, p, logger_ta, logger_rh FROM meteoseries WHERE stationID='STATION_NUMBER' AND TimeStamp>='START_DATE' AND TimeStamp<='END_DATE' ORDER BY timestamp ASC
+	std::string Query2( MySQLQueryMeteoData );
+	const std::string str1( "STATION_NAME" );
+	const std::string str2( "STATION_NUMBER" );
+	const std::string str3( "START_DATE" );
+	const std::string str4( "END_DATE" );
+	Query2.replace( Query2.find(str1), str1.length(), stat_abk ); // set 1st variable's value
+	Query2.replace( Query2.find(str2), str2.length(), stao_nr ); // set 2nd variable's value
+	Query2.replace( Query2.find(str3), str3.length(), sDate ); // set 3rd variable's value
+	Query2.replace( Query2.find(str4), str4.length(), eDate ); // set 4th variable's value
+
+	MYSQL *conn2 = mysql_init(nullptr);
+	if (!mysql_real_connect(conn2, mysqlhost.c_str(), mysqluser.c_str(), mysqlpass.c_str(), mysqldb.c_str(), 0, nullptr, 0)) {
+		throw IOException("Could not initiate connection to Mysql server "+mysqlhost, AT);
+	}
+
+	if (mysql_query(conn2, Query2.c_str())) {
+		throw IOException("Query2 \""+Query2+"\" failed to execute", AT);
+	}
+
+	MYSQL_RES *res2 = mysql_use_result(conn2);
+	const unsigned int column_no2 = mysql_num_fields(res2);
+	MYSQL_ROW row2;
+	while ( ( row2= mysql_fetch_row(res2) ) != nullptr ) {
+		std::vector<std::string> vecData;
+		for (unsigned int ii=0; ii<column_no2; ii++) {
+			std::string row_02;
+				if (!row2[ii]){
+					IOUtils::convertString(row_02,"-999.0");// HARD CODED :(
+				}else{
+					IOUtils::convertString(row_02, row2[ii]);
+				 }
+			vecData.push_back(row_02);
+		}
+		if (fullStation) {
+			for (unsigned int ii=0; ii<static_cast<unsigned int>(vecHTS1.size()); ii++) {
+				vecData.push_back(vecHTS1.at(ii));
+			}
+		}
+		vecMeteoData.push_back(vecData);
+	}
+
+	mysql_free_result(res2);
+	mysql_close(conn2);
+	return fullStation;
+}
+
+} //namespace
diff --git a/Source/meteoio/meteoio/plugins/WWCSIO.h b/Source/meteoio/meteoio/plugins/WWCSIO.h
new file mode 100644
index 00000000..d1a83271
--- /dev/null
+++ b/Source/meteoio/meteoio/plugins/WWCSIO.h
@@ -0,0 +1,75 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+/***********************************************************************************/
+/*  Copyright 2014 Snow and Avalanche Study Establishment    SASE-CHANDIGARH       */
+/***********************************************************************************/
+/* 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/>.
+*/
+#ifndef WWCSIO_H
+#define WWCSIO_H
+
+#include <meteoio/IOInterface.h>
+
+#include <string>
+
+namespace mio {
+
+/**
+ * @class WWCSIO_H
+ * @brief This is the plugin required to get meteorological data from the WWCS database.
+ *
+ * @ingroup plugins
+ * @author Mathias Bavay
+ * @date   2022-02-24
+ */
+class WWCSIO : public IOInterface {
+	public:
+		WWCSIO(const std::string& configfile);
+		WWCSIO(const WWCSIO&);
+		WWCSIO(const Config& cfgreader);
+
+		virtual void readStationData(const Date& date, std::vector<StationData>& vecStation);
+		virtual void readMeteoData(const Date& dateStart, const Date& dateEnd,
+		                           std::vector< std::vector<MeteoData> >& vecMeteo);
+
+	private:
+		void readConfig();
+		void readStationIDs(std::vector<std::string>& vecStationID) const;
+		static void parseStationID(const std::string& stationID, std::string& stnAbbrev, std::string& stnNumber);
+		void getStationMetaData(const std::string& stat_abk, const std::string& stao_nr,const std::string& sqlQuery,
+		                        std::vector<std::string>& vecMetaData);
+		void readStationMetaData();
+		void readData(const Date& dateStart, const Date& dateEnd, std::vector< std::vector<MeteoData> >& vecMeteo,
+		              const size_t& stationindex, const std::vector<StationData>& vecMeta) const;
+		static void convertUnits(MeteoData& meteo);
+		void parseDataSet(const std::vector<std::string>& i_meteo, MeteoData& md) const;
+		bool getStationData(const std::string& stat_abk, const std::string& stao_nr,const Date& dateS,
+		                    const Date& dateE,const std::vector<std::string>& vecHTS1,
+		                    std::vector< std::vector<std::string> >& vecMeteoData) const;
+
+		const Config cfg;
+		std::vector<std::string> vecStationIDs;
+		std::vector<StationData> vecStationMetaData;
+		std::string mysqlhost, mysqldb, mysqluser, mysqlpass;
+		std::string coordin, coordinparam, coordout, coordoutparam; //projection parameters
+		double in_dflt_TZ, out_dflt_TZ;
+
+		static const double plugin_nodata; //plugin specific nodata value, e.g. -999
+		static const std::string MySQLQueryStationMetaData;
+		static const std::string MySQLQueryMeteoData;
+};
+
+} //namespace
+#endif
+
diff --git a/Source/snowpack/snowpack/DataClasses.cc b/Source/snowpack/snowpack/DataClasses.cc
index ad563849..f94db0c1 100644
--- a/Source/snowpack/snowpack/DataClasses.cc
+++ b/Source/snowpack/snowpack/DataClasses.cc
@@ -370,12 +370,12 @@ void SurfaceFluxes::collectSurfaceFluxes(const BoundCond& Bdata,
 	// 2) Long wave fluxes.
 	lw_out += Bdata.lw_out;
 	lw_net += Bdata.lw_net;
-	if (Mdata.lw_net != IOUtils::nodata) {
+	if (Mdata.lw_net == IOUtils::nodata) {
 		// Default
 		lw_in  +=  Atmosphere::blkBody_Radiation(Mdata.ea, Mdata.ta);
 	} else {
 		// NET_LW provided
-		lw_in  +=  lw_net + lw_out;
+		lw_in  +=  Bdata.lw_net + Bdata.lw_out;
 	}
 
 	// 3) Turbulent fluxes.