Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/SnowpackInterfaceWorker.cc

1/***********************************************************************************/
2/* Copyright 2009-2015 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3/***********************************************************************************/
4/* This file is part of Alpine3D.
5 Alpine3D is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Alpine3D is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with Alpine3D. If not, see <http://www.gnu.org/licenses/>.
17*/
18#include <alpine3d/SnowpackInterfaceWorker.h>
19#include <alpine3d/AlpineMain.h>
20
21using namespace std;
22using namespace mio;
23
24//This is a function to retrieve soil temperature at a given depth (measured from
25//the soil surface, not from the snow surface) for output purposes.
26inline double getSoilTemperature(const SnowStation& pixel, const double& depth)
27{
28if(pixel.getNumberOfNodes() == 0) {
29return IOUtils::nodata;
30} else if(depth < 1e-5) {//we want the temperature at the soil surface
31return pixel.Ndata[pixel.SoilNode].T;
32} else if(depth >= pixel.Ground) {//we want the temperature below the lowest node
33return pixel.Ndata.front().T;
34} else {
35const double height( pixel.Ground - depth );
36
37//Looking for the node located directly above the specified depth
38size_t iNode(0);
39while (pixel.Ndata[iNode].z < height)
40++iNode;
41
42const NodeData& aboveNode( pixel.Ndata[iNode] );
43if(aboveNode.z == height) {//exact match
44return aboveNode.T;
45} else {//compute linear interpolation between the two nodes
46const NodeData& belowNode( pixel.Ndata[iNode-1] );
47const double a = (aboveNode.T - belowNode.T)/(aboveNode.z - belowNode.z);
48return a*(height - belowNode.z) + belowNode.T;
49}
50}
51}
52
53//This is a function to retrieve a given parameter at a given depth in the snow
54//for output purposes.
55inline double getValueAtDepth(const SnowStation& pixel, const double& depth)
56{
57const size_t nNodes = pixel.getNumberOfNodes();
58if (nNodes==0) return mio::IOUtils::nodata; //no nodes -> nodata
59
60const size_t max_ii = nNodes - 1;
61const double offset = pixel.Ground;
62const double snow_height = pixel.Ndata[max_ii].z - offset;
63const double height = snow_height - depth;
64
65if (height<0.) return mio::IOUtils::nodata; //not deep enough -> nodata
66if (depth==0.) return pixel.Ndata[max_ii].T; //if depth==0 -> surface
67
68//looking for the first node less than depth.
69//At this point, we know that we have >= than depth snow
70size_t ii = 0;
71while ( (pixel.Ndata[ii].z-offset) <= height ) {
72ii++;
73}
74
75if ( (pixel.Ndata[ii].z-offset) == height) return pixel.Ndata[ii].T; //we found it exactly
76
77//compute linear interpolation between the two points y = az+b
78const double a = (pixel.Ndata[ii].T-pixel.Ndata[ii-1].T) / (pixel.Ndata[ii].z-pixel.Ndata[ii-1].z);
79const double b = pixel.Ndata[ii-1].T;
80const double z = height - (pixel.Ndata[ii-1].z-offset);
81return a*z + b;
82}
83
84// Returns the average snow temperature in the top x cm, or, when snow depth is less than x cm, the average snow temperature.
85inline double getAvgAtDepth(const SnowStation& pixel, const double& depth, const bool& is_temp)
86{
87
88double thickness = 0.;
89double sum = 0.;
90if (pixel.getNumberOfElements() > pixel.SoilNode) {
91size_t i = pixel.getNumberOfElements();
92while (i-- > pixel.SoilNode) {
93const double val = (is_temp)? pixel.Edata[i].Te : pixel.Edata[i].Rho;
94if (thickness < depth) {
95sum += pixel.Edata[i].L * val;
96thickness += pixel.Edata[i].L;
97} else { //we found the first element that is too deep
98if (thickness!=depth) { //just in case, we might have found an element exactly at "depth"
99const double thick_contrib = depth - thickness;
100sum += thick_contrib * val;
101thickness += thick_contrib;
102}
103break;
104}
105}
106}
107
108const double value = (thickness == 0.) ? (IOUtils::nodata) : (sum / thickness);
109return value;
110}
111
112/******************************************************************************
113 * Constructor / Destructor
114 ******************************************************************************/
115
116/**
117 * @brief Constructs and initialise Snowpack Interface Worker. Create one worker
118 * and init values for slice which the worker need to simulate it slice.
119 * @param io_cfg configuration to pass to Snowpack
120 * @param dem_in gives the demographic Data. Also tetermines size and position of the slice
121 * @param landuse_in gives the landuse Data. Also tetermines size and position of the landuse for slice
122 * @param pts_in gives the spezial points. For this points more output is done then for the others. Calcualtion is the same.
123 * @param snow_stations gives a vector of pointers to the SnowStation objects relevant for this thread
124 * @param snow_stations_coord provide the (ii,jj) coordinates of each snow station. This is necessary because the slices might be irregular
125 * @param offset_in gives the offsett on X for this slice (needed to read data and error messages)
126 */
127SnowpackInterfaceWorker::SnowpackInterfaceWorker(const mio::Config& io_cfg,
128 const mio::DEMObject& dem_in,
129 const mio::Grid2DObject& landuse_in,
130 const std::vector< std::pair<size_t,size_t> >& pts_in,
131 const std::vector<SnowStation*>& snow_stations,
132 const std::vector<std::pair<size_t,size_t> >& snow_stations_coord,
133 const size_t offset_in)
134 : sn_cfg(io_cfg), sn(sn_cfg), meteo(sn_cfg), stability(sn_cfg, false), sn_techsnow(sn_cfg), dem(dem_in),
135 dimx(dem.getNx()), dimy(dem.getNy()), offset(offset_in), SnowStations(snow_stations), SnowStationsCoord(snow_stations_coord),
136 isSpecialPoint(snow_stations.size(), false), landuse(landuse_in), store(dem_in, 0.), grids(), snow_pixel(), meteo_pixel(),
137 surface_flux(), soil_temp_depths(), calculation_step_length(0.), height_of_wind_value(0.),
138 snow_temp_depth(IOUtils::nodata), snow_avg_temp_depth(IOUtils::nodata), snow_avg_rho_depth(IOUtils::nodata),
139 useDrift(false), useEBalance(false), useCanopy(false)
140{
141
142sn_cfg.getValue("CALCULATION_STEP_LENGTH", "Snowpack", calculation_step_length);
143sn_cfg.getValue("HEIGHT_OF_WIND_VALUE", "Snowpack", height_of_wind_value); //currently unused
144sn_cfg.getValue("CANOPY", "Snowpack", useCanopy);
145io_cfg.getValue("SNOW_TEMPERATURE_DEPTH", "Output", snow_temp_depth, IOUtils::nothrow);
146io_cfg.getValue("SNOW_AVG_TEMPERATURE_DEPTH", "Output", snow_avg_temp_depth, IOUtils::nothrow);
147io_cfg.getValue("SNOW_AVG_DENSITY_DEPTH", "Output", snow_avg_rho_depth, IOUtils::nothrow);
148
149//create the vector of output grids
150std::vector<std::string> params = sn_cfg.get("GRIDS_PARAMETERS", "output");
151if (snow_temp_depth!=IOUtils::nodata) params.push_back("TSNOW");
152if (snow_avg_temp_depth!=IOUtils::nodata) params.push_back("TSNOW_AVG");
153if (snow_avg_rho_depth!=IOUtils::nodata) params.push_back("RHOSNOW_AVG");
154
155//handle the soil temperatures
156io_cfg.getValue("SOIL_TEMPERATURE_DEPTHS", "Output", soil_temp_depths, IOUtils::nothrow);
157const unsigned short max_Tsoil( SnGrids::lastparam - SnGrids::TSOIL1 + 1 );
158if (soil_temp_depths.size()>max_Tsoil)
159throw InvalidArgumentException("Too many soil temperatures requested", AT);
160for (size_t ii=0; ii<soil_temp_depths.size(); ii++) {
161params.push_back( "TSOIL"+mio::IOUtils::toString(ii+1) );
162}
163uniqueOutputGrids(params);
164initGrids(params);
165
166const CurrentMeteo meteoPixel; //this is only necessary in order to have something for fillGrids()
167const SurfaceFluxes surfaceFlux; //this is only necessary in order to have something for fillGrids()
168
169for(size_t ii = 0; ii < SnowStationsCoord.size(); ++ii) {
170size_t ix = SnowStationsCoord.at(ii).first;
171size_t iy = SnowStationsCoord.at(ii).second;
172if (SnowpackInterfaceWorker::skipThisCell(landuse(ix,iy), dem(ix,iy))) { //skip nodata cells as well as water bodies, etc
173if (!pts_in.empty() && is_special(pts_in, ix, iy)){
174std::cout << "[W] POI (" << ix << "," << iy << ") will be skipped (nodata, water body, etc)" << std::endl;
175}
176continue;
177}
178
179if (SnowStations[ii]==NULL) continue; //for safety: skip cells initialized with NULL
180const SnowStation &snowPixel = *SnowStations[ii];
181fillGrids(ix, iy, meteoPixel, snowPixel, surfaceFlux);
182if (!pts_in.empty()) { //is it a special point?
183isSpecialPoint[ii] = is_special(pts_in, ix, iy);
184}
185}
186}
187
188SnowpackInterfaceWorker::~SnowpackInterfaceWorker()
189{
190//sorry for this cryptic syntax, this is to guarantee execution order with the "," operator
191while (!SnowStations.empty())
192((SnowStations.back()!=NULL)? delete SnowStations.back() : (void)0) , SnowStations.pop_back();
193}
194
195/******************************************************************************
196 * Methods that have to do with output
197 ******************************************************************************/
198
199/** @brief Make sure all requested grids only appear once
200 * @param output_grids vector of requeste grids to sort and filter
201 */
202void SnowpackInterfaceWorker::uniqueOutputGrids(std::vector<std::string>& output_grids)
203{
204for (size_t ii = 0; ii<output_grids.size(); ++ii)
205IOUtils::toUpper(output_grids[ii]);
206
207std::sort(output_grids.begin(), output_grids.end());
208const std::vector<std::string>::iterator it = std::unique(output_grids.begin(), output_grids.end());
209output_grids.resize(std::distance(output_grids.begin(),it));
210}
211
212/** @brief Initialize and add to the grid map the requested grids
213 * @param params string representation of the grids to add
214 */
215void SnowpackInterfaceWorker::initGrids(std::vector<std::string>& params)
216{
217for (size_t ii = 0; ii<params.size(); ++ii) {
218IOUtils::toUpper(params[ii]); //make sure all parameters are upper case
219
220const size_t param_idx = SnGrids::getParameterIndex( params[ii] );
221if (param_idx==IOUtils::npos)
222throw UnknownValueException("Unknow meteo grid '"+params[ii]+"' selected for gridded output", AT);
223
224const std::map< SnGrids::Parameters, mio::Grid2DObject >::const_iterator it = grids.find( static_cast<SnGrids::Parameters>(param_idx) );
225if (it==grids.end()) { //the parameter did not already exist, adding it
226Grid2DObject tmp_grid(dem, IOUtils::nodata);
227grids[ static_cast<SnGrids::Parameters>(param_idx) ] = tmp_grid;
228}
229}
230}
231
232void SnowpackInterfaceWorker::getOutputSpecialPoints(std::vector<SnowStation*>& ptr_snow_pixel,
233 std::vector<CurrentMeteo*>& ptr_meteo_pixel,
234 std::vector<SurfaceFluxes*>& ptr_surface_flux)
235{
236for (size_t ii=0; ii<snow_pixel.size(); ii++) {
237ptr_snow_pixel.push_back( &snow_pixel[ii] );
238ptr_meteo_pixel.push_back( &meteo_pixel[ii] );
239ptr_surface_flux.push_back( &surface_flux[ii] );
240}
241}
242
243void SnowpackInterfaceWorker::clearSpecialPointsData()
244{
245snow_pixel.clear();
246meteo_pixel.clear();
247surface_flux.clear();
248}
249
250/**
251 * @brief method that returns SnowCover files for specific date.
252 * @param snow_station vector that will be filled with the SnowStation data for each pixels
253 */
254void SnowpackInterfaceWorker::getOutputSNO(std::vector<SnowStation*>& snow_station) const
255{
256for(size_t i = 0; i < SnowStationsCoord.size(); ++i) {
257size_t ix = SnowStationsCoord.at(i).first;
258size_t iy = SnowStationsCoord.at(i).second;
259if (SnowpackInterfaceWorker::skipThisCell(landuse(ix,iy), dem(ix,iy))) continue; //skip nodata cells as well as water bodies, etc
260const size_t index_SnowStation = i;
261if (SnowStations[index_SnowStation]==NULL) continue; //for safety: skipped cells were initialized with NULL
262#pragma omp critical (snow_station_lock)
263snow_station.push_back( SnowStations[index_SnowStation] );
264}
265}
266
267/****************************************************************************
268 * Getters and Setters
269 ****************************************************************************/
270
271/**
272 * @brief Method that the Master can search the neded data (in grids) from Worker (Pull from client)
273 * @param param says which grid param the Master wants to have
274 * @return the 2D output grid, which gives back the data to the master
275 */
276mio::Grid2DObject SnowpackInterfaceWorker::getGrid(const SnGrids::Parameters& param) const
277{
278const std::map< SnGrids::Parameters, mio::Grid2DObject >::const_iterator it = grids.find(param);
279if (it==grids.end()) {
280throw UnknownValueException("Undeclared grid '"+SnGrids::getParameterName(param)+"' requested from SnowpackInterfaceWorker", AT);
281}
282
283return it->second;
284}
285
286/**
287 * @brief Retrieve one point (ii,jj) from the specified grid.
288 * @param param says which grid param the Master wants to have
289 * @param ii ii index
290 * @param jj jj index
291 * @return grid value at point (ii,jj) for parameter param
292 */
293double& SnowpackInterfaceWorker::getGridPoint(const SnGrids::Parameters& param, const size_t& ii, const size_t& jj)
294{
295const std::map< SnGrids::Parameters, mio::Grid2DObject >::const_iterator it( grids.find(param) );
296if (it==grids.end())
297throw UnknownValueException("Undeclared grid '"+SnGrids::getParameterName(param)+"' requested from SnowpackInterfaceWorker", AT);
298
299return grids[ param ](ii,jj);
300}
301
302/**
303 * @brief Fill all the grids stored in the **grids** map of 2D grids with all the necessary values for the provided pixel.
304 * Before calling this method, make sure that snowPixel is not NULL (as it could be for pixels that should be skipped)
305 * @param ii ii index
306 * @param jj jj index
307 * @param meteoPixel meteorological forcing for the pixel
308 * @param snowPixel canopy/snow/soil stratigraphy information
309 * @param surfaceFlux surface fluxes for the pixel
310 */
311void SnowpackInterfaceWorker::fillGrids(const size_t& ii, const size_t& jj, const CurrentMeteo& meteoPixel, const SnowStation& snowPixel, const SurfaceFluxes& surfaceFlux)
312{
313std::map< SnGrids::Parameters, mio::Grid2DObject >::iterator it;
314for (it=grids.begin(); it!=grids.end(); ++it) {
315double value = IOUtils::nodata;
316switch (it->first) {
317case SnGrids::TA:
318value = meteoPixel.ta; break;
319case SnGrids::RH:
320value = meteoPixel.rh; break;
321case SnGrids::VW:
322value = meteoPixel.vw; break;
323case SnGrids::DW:
324value = meteoPixel.dw; break;
325case SnGrids::ISWR:
326value = meteoPixel.iswr; break;
327case SnGrids::ISWR_DIFF:
328value = (useEBalance)? meteoPixel.diff : IOUtils::nodata; break;
329case SnGrids::ISWR_DIR:
330value = (useEBalance)? meteoPixel.iswr - meteoPixel.diff : IOUtils::nodata; break;
331case SnGrids::ILWR:
332value = Atmosphere::blkBody_Radiation(meteoPixel.ea, meteoPixel.ta); break;
333case SnGrids::HS:
334value = (snowPixel.cH - snowPixel.Ground) / snowPixel.cos_sl; break; //slope2horiz
335case SnGrids::PSUM:
336value = meteoPixel.psum; break;
337case SnGrids::PSUM_PH:
338value = meteoPixel.psum_ph; break;
339case SnGrids::PSUM_TECH:
340value = meteoPixel.psum_tech; break;
341case SnGrids::TSS:
342if (!useCanopy || snowPixel.Cdata.zdispl < 0.) {
343value = snowPixel.Ndata.back().T;
344} else { // radiative surface temperature above canopy
345value = pow(snowPixel.Cdata.rlwrac/Constants::stefan_boltzmann, 0.25);
346}
347break;
348case SnGrids::TSG:
349value = snowPixel.Ndata[snowPixel.SoilNode].T; break;
350case SnGrids::TS0:
351value = (snowPixel.SoilNode>0)? snowPixel.Ndata[snowPixel.SoilNode].T : IOUtils::nodata; break;
352case SnGrids::TSNOW:
353value = getValueAtDepth(snowPixel, snow_temp_depth); break;
354case SnGrids::TSNOW_AVG:
355value = getAvgAtDepth(snowPixel, snow_temp_depth, true); break;
356case SnGrids::RHOSNOW_AVG:
357value = getAvgAtDepth(snowPixel, snow_temp_depth, false); break;
358case SnGrids::SWE:
359value = surfaceFlux.mass[SurfaceFluxes::MS_TOTALMASS] / snowPixel.cos_sl; break; //slope2horiz
360case SnGrids::RSNO: {
361const double hs = (snowPixel.cH - snowPixel.Ground);
362value = (hs>0.)? snowPixel.mass_sum / hs : IOUtils::nodata;
363break; }
364case SnGrids::TOP_ALB:
365if (!useCanopy || snowPixel.Cdata.zdispl < 0.) {
366value = snowPixel.Albedo;
367} else { // above canopy albedo
368value = snowPixel.Cdata.totalalb;
369}
370break;
371case SnGrids::SURF_ALB:
372value = snowPixel.Albedo; break;
373case SnGrids::SP:
374value = (!snowPixel.Edata.empty())? snowPixel.Edata.back().sp : IOUtils::nodata; break;
375case SnGrids::RB:
376value = (!snowPixel.Edata.empty())? snowPixel.Edata.back().rb : IOUtils::nodata; break;
377case SnGrids::RG:
378value = (!snowPixel.Edata.empty())? snowPixel.Edata.back().rg : IOUtils::nodata; break;
379case SnGrids::N3:
380value = (!snowPixel.Edata.empty())? snowPixel.Edata.back().N3 : IOUtils::nodata; break;
381case SnGrids::MS_SNOWPACK_RUNOFF:
382value = surfaceFlux.mass[SurfaceFluxes::MS_SNOWPACK_RUNOFF] / snowPixel.cos_sl; break;
383case SnGrids::MS_SOIL_RUNOFF:
384value = surfaceFlux.mass[SurfaceFluxes::MS_SOIL_RUNOFF] / snowPixel.cos_sl; break;
385case SnGrids::MS_WATER:
386value = surfaceFlux.mass[SurfaceFluxes::MS_WATER] / snowPixel.cos_sl; break;
387case SnGrids::SFC_SUBL:
388value = -surfaceFlux.mass[SurfaceFluxes::MS_SUBLIMATION] / snowPixel.cos_sl; break; //slope2horiz
389case SnGrids::STORE:
390value = store(ii,jj); break;
391case SnGrids::GLACIER:
392value = (!snowPixel.isGlacier(true))? 1. : IOUtils::nodata; break; //glaciated pixels receive IOUtils::nodata
393case SnGrids::GLACIER_EXPOSED:
394value = (!snowPixel.isGlacier(false))? 1. : IOUtils::nodata; break; //glaciated pixels receive IOUtils::nodata
395case SnGrids::ET:
396value = -(surfaceFlux.mass[SurfaceFluxes::MS_SUBLIMATION]+surfaceFlux.mass[SurfaceFluxes::MS_EVAPORATION])/snowPixel.cos_sl; //slope2horiz
397// Add part from Canopy
398value += useCanopy?(snowPixel.Cdata.transp+snowPixel.Cdata.intevap)/snowPixel.cos_sl:0; //slope2horiz
399break;
400default:
401if (it->first>=SnGrids::TSOIL1 && it->first<=SnGrids::lastparam) //dealing with soil temperatures
402{
403value = (soil_temp_depths.empty())? IOUtils::nodata : getSoilTemperature(snowPixel, soil_temp_depths[ it->first - SnGrids::TSOIL1 ]);
404}
405else
406{
407throw InvalidArgumentException("Invalid parameter requested", AT);
408}
409}
410it->second(ii,jj) = value;
411}
412}
413
414/**
415 * @brief method which prepares all data for simulation and then access correctly the
416 * Snowpack model interfaces.
417 * @param date current simulation time step
418 * @param psum precipitation grid (kg m-2)
419 * @param psum_ph precipitation phase grid (between 0 and 1)
420 * @param psum_tech technical precipitation grid (kg m-2)
421 * @param rh relative humidity grid (% or 1)
422 * @param ta air temperature grid (K)
423 * @param vw wind velocity grid (m s-1)
424 * @param dw wind direction grid (degrees north)
425 * @param mns map of the Precipitation (mm/h) HACK get this map only if per pull from Master if Drift is used !!
426 * @param shortwave incoming shortwave radiation grid (W m-2)
427 * @param diffuse diffuse radiation from the sky grid (W m-2)
428 * @param longwave incoming longwave grid (W m-2)
429 * @param solarElevation solar elevation (in dec)
430 */
431void SnowpackInterfaceWorker::runModel(const mio::Date &date,
432 const mio::Grid2DObject &psum,
433 const mio::Grid2DObject &psum_ph,
434 const mio::Grid2DObject &psum_tech,
435 const mio::Grid2DObject &rh,
436 const mio::Grid2DObject &ta,
437 const mio::Grid2DObject &vw,
438 const mio::Grid2DObject &dw,
439 const mio::Grid2DObject &mns,
440 const mio::Grid2DObject &shortwave,
441 const mio::Grid2DObject &diffuse,
442 const mio::Grid2DObject &longwave,
443 const double solarElevation)
444{
445const Meteo::ATM_STABILITY USER_STABILITY = meteo.getStability();
446const std::string bcu_watertransportmodel_snow = sn_cfg.get("WATERTRANSPORTMODEL_SNOW", "SnowpackAdvanced");
447const std::string bcu_watertransportmodel_soil = sn_cfg.get("WATERTRANSPORTMODEL_SOIL", "SnowpackAdvanced");
448const std::string bcu_reduce_n_elements = sn_cfg.get("REDUCE_N_ELEMENTS", "SnowpackAdvanced");
449const std::string bcu_adjust_height_of_meteo= sn_cfg.get("ADJUST_HEIGHT_OF_METEO_VALUES", "SnowpackAdvanced");
450const std::string bcu_adjust_height_of_wind = sn_cfg.get("ADJUST_HEIGHT_OF_WIND_VALUE", "SnowpackAdvanced");
451
452CurrentMeteo meteoPixel(sn_cfg);
453meteoPixel.date = date;
454meteoPixel.elev = solarElevation*Cst::to_rad; //HACK: Snowpack uses RAD !!!!!
455
456// make SN calculations....
457for(size_t i = 0; i < SnowStationsCoord.size(); ++i) {
458size_t ix = SnowStationsCoord.at(i).first;
459size_t iy = SnowStationsCoord.at(i).second;
460if (SnowpackInterfaceWorker::skipThisCell(landuse(ix,iy), dem(ix,iy))) continue; //skip nodata cells as well as water bodies, etc
461const size_t index_SnowStation = i;
462if (SnowStations[index_SnowStation]==NULL) continue; //for safety: skipped cells were initialized with NULL
463SnowStation &snowPixel = *SnowStations[index_SnowStation];
464const bool isGlacier = snowPixel.isGlacier(false);
465
466//In case of ice and firn pixels, use BUCKET model for water transport:
467const int land = (round_landuse(landuse(ix,iy)) - 10000) / 100;
468if (land==13 || land==14) {
469sn_cfg.addKey("WATERTRANSPORTMODEL_SNOW", "SnowpackAdvanced", "BUCKET");
470sn_cfg.addKey("WATERTRANSPORTMODEL_SOIL", "SnowpackAdvanced", "BUCKET");
471}
472// In case of glacier pixel, remove meteo height correction and try to merge elemnts
473if (land==14) {
474sn_cfg.addKey("REDUCE_N_ELEMENTS", "SnowpackAdvanced", "TRUE");
475sn_cfg.addKey("ADJUST_HEIGHT_OF_METEO_VALUES","SnowpackAdvanced","FALSE");
476sn_cfg.addKey("ADJUST_HEIGHT_OF_WIND_VALUE","SnowpackAdvanced","FALSE");
477}
478// Set curent meteo variables from 2D fields to single pixel
479const double previous_albedo = getGridPoint(SnGrids::TOP_ALB, ix, iy);
480meteoPixel.rh = rh(ix,iy);
481meteoPixel.ta = ta(ix,iy);
482meteoPixel.vw = vw(ix,iy);
483meteoPixel.dw = dw(ix,iy);
484meteoPixel.iswr = shortwave(ix,iy);
485meteoPixel.rswr = previous_albedo*meteoPixel.iswr;
486meteoPixel.tss = snowPixel.Ndata[snowPixel.getNumberOfElements()].T; //we use previous timestep value
487meteoPixel.ts0 = (snowPixel.SoilNode>0)? snowPixel.Ndata[snowPixel.SoilNode].T : meteoPixel.tss; //we use previous timestep value
488meteoPixel.ea = Atmosphere::blkBody_Emissivity(longwave(ix,iy), meteoPixel.ta); //to be consistent with Snowpack
489meteoPixel.psum_ph = psum_ph(ix,iy);
490meteoPixel.psum_tech = psum_tech(ix, iy);
491meteoPixel.hs = IOUtils::nodata;
492if ( (meteoPixel.tss != IOUtils::nodata && meteoPixel.tss<=100) || (meteoPixel.ts0 != IOUtils::nodata && meteoPixel.ts0<=100) ) {
493cout << "[E] Pixel (" << ix << "," << iy << ") too cold! tss=" << meteoPixel.tss << " ts0=" << meteoPixel.ts0 << std::endl;
494}
495
496// Now determine perpendicular to slope mass balance from drift or precip
497// Note that psum must be in mm/h and SNOWPACK will converted it to mm/CALCULTION_STEP_LENGTH
498if ( useDrift && mns(ix,iy)!=IOUtils::nodata ) {
499double drift_mass = mns(ix,iy);
500// For extreme terrain, the drift might go crazy: limit deposition / erosion
501if (fabs(drift_mass)>37.) {
502cout << "crazy drift at (" << ix << "," << iy << ") = " << drift_mass << " mm/h\n";
503drift_mass = (drift_mass>0.)? 37. : -37.; //set to deposition limit
504}
505
506meteoPixel.psum = drift_mass;
507} else {
508meteoPixel.psum= psum(ix,iy) * snowPixel.cos_sl; //horiz2slope
509}
510
511if (useEBalance) meteoPixel.diff = diffuse(ix,iy);
512// Reset Canopy Surface Data to 0 before calling Meteo and Canopy module
513if (useCanopy) snowPixel.Cdata.initializeSurfaceExchangeData();
514//if the current pixel contains soil layers, then activate soil modeling in Snowpack
515sn.setUseSoilLayers( snowPixel.hasSoilLayers() ); //TODO: when Snowpack runs, it should check it
516
517// exposed glacier special case
518if (isGlacier) {
519const std::string tmp_sw_mode = sn_cfg.get("SW_MODE", "Snowpack");
520if (tmp_sw_mode == "BOTH") {
521//switch to glacier albedo (when sw_mode != BOTH, the calculation internally relies on SnLaws and should not be overwritten here!)
522snowPixel.Albedo = Constants::glacier_albedo;
523}
524if (meteoPixel.ta>IOUtils::C_TO_K(5.)) {
525//switch to STABLE atmosphere on glacier if TA>5°C
526meteo.setStability(Meteo::MO_HOLTSLAG);
527}
528}
529
530//compute ustar, psi_s, z0
531meteo.compMeteo(meteoPixel, snowPixel, true);
532SurfaceFluxes surfaceFlux;
533// run snowpack model itself
534double dIntEnergy = 0.; //accumulate the dIntEnergy over the snowsteps
535const unsigned int nr_snowsteps = (unsigned int)(dt_main/M_TO_S(calculation_step_length));
536for (unsigned int snowsteps = 0; snowsteps < nr_snowsteps; snowsteps++) {
537/* Update the store variable */
538/* david: why += ? isnt store reset every timestep ? */
539/* Michi: This is exactly the point, store is only set to zero when enough precipitation
540has accumulated to create at least one new (finite element) layer */
541
542if (snowsteps == 0) meteoPixel.psum /= nr_snowsteps;
543store(ix,iy) += meteoPixel.psum;
544
545try {
546BoundCond Bdata;
547sn.runSnowpackModel(meteoPixel, snowPixel, store(ix,iy), Bdata, surfaceFlux);
548surfaceFlux.collectSurfaceFluxes(Bdata, snowPixel, meteoPixel);
549dIntEnergy += snowPixel.dIntEnergy; //it is reset at every new call to runSnowpackModel
550meteoPixel.hs = snowPixel.cH - snowPixel.Ground; //do not reproject here, otherwise Snowpack outputs would get messed up
551} catch (const std::bad_alloc&) { //don't try anything fancy when running low on memory
552const int lus =(int)floor( landuse(ix,iy));
553const double slope2horiz = (1. / snowPixel.cos_sl);
554const double snow = (snowPixel.cH - snowPixel.Ground) * slope2horiz;
555cout << "[E] Could not allocate memory in Snowpack for cell (" << ix << "," << iy << ") LUS=" << lus << " ";
556cout << "with " << std::fixed << std::setprecision(4) << snow << " m snow in " << snowPixel.getNumberOfElements()-snowPixel.SoilNode << "/" << snowPixel.getNumberOfElements() << " elements.\n";
557fflush( stdout );
558throw;
559} catch (const std::exception& e) {
560const int lus =(int)floor( landuse(ix,iy));
561const double slope2horiz = (1. / snowPixel.cos_sl);
562const double snow = (snowPixel.cH - snowPixel.Ground) * slope2horiz;
563if (nr_snowsteps > 1) {
564surfaceFlux.multiplyFluxes(1./nr_snowsteps);
565if (useCanopy)
566snowPixel.Cdata.multiplyFluxes(1./nr_snowsteps);
567meteoPixel.psum *= nr_snowsteps;
568snowPixel.dIntEnergy = dIntEnergy;
569}
570
571ostringstream ss;
572ss << "[E] Snowpack exception: " << e.what() << "\n";
573ss << "[E] at cell (" << ix << "," << iy << ") LUS=" << lus << " ";
574ss << "with " << std::fixed << std::setprecision(4) << snow << " m snow in " << snowPixel.getNumberOfElements()-snowPixel.SoilNode << "/" << snowPixel.getNumberOfElements() << " elements.\n";
575if (isSpecialPoint[index_SnowStation])
576gatherSpecialPoints(meteoPixel, snowPixel, surfaceFlux); //gather special point data, in order to get as much information as possible
577throw IOException(ss.str(), AT);
578}
579
580// Fill the surfaceFlux.mass variable for output
581surfaceFlux.mass[SurfaceFluxes::MS_HNW] += meteoPixel.psum;
582}
583
584//some variables are now wrong if we ran multiple Snowpack steps -> recompute them!
585if (nr_snowsteps > 1) {
586surfaceFlux.multiplyFluxes(1./nr_snowsteps);
587if (useCanopy)
588snowPixel.Cdata.multiplyFluxes(1./nr_snowsteps);
589meteoPixel.psum *= nr_snowsteps;
590snowPixel.dIntEnergy = dIntEnergy;
591}
592
593//switch stability back to normal if it was changed
594if (meteo.getStability()!=USER_STABILITY) meteo.setStability(USER_STABILITY);
595//if the glacier is still exposed, force the albedo back to glacier albedo
596if (isGlacier) {
597const std::string tmp_sw_mode = sn_cfg.get("SW_MODE", "Snowpack");
598if (tmp_sw_mode == "BOTH") {
599//switch to glacier albedo (when sw_mode != BOTH, the calculation internally relies on SnLaws and should not be overwritten here!)
600surfaceFlux.pAlbedo = snowPixel.Albedo = Constants::glacier_albedo;
601}
602}
603if (!std::isfinite( getGridPoint(SnGrids::TOP_ALB, ix, iy) )) {
604//if the albedo is nan, infinity, etc reset it to its previous
605//value to try to rescue the pixel...
606cerr << "[E] pixel (" << ix << "," << iy << ") found with a nan/infinit albedo ["<< getGridPoint(SnGrids::TOP_ALB, ix, iy) <<"]; reseting to " << previous_albedo << std::endl;
607getGridPoint(SnGrids::TOP_ALB, ix, iy) = previous_albedo;
608}
609
610try{
611stability.checkStability(meteoPixel, snowPixel);
612} catch(...) {
613snowPixel.S_4 = Stability::max_stability; //nothing else to do...
614}
615
616surfaceFlux.mass[SurfaceFluxes::MS_TOTALMASS] = 0.0;
617const std::vector<ElementData> EMS( snowPixel.Edata );
618for (size_t e=0; e<EMS.size(); e++) {
619if (EMS[e].theta[SOIL] <= 0.) {
620surfaceFlux.mass[SurfaceFluxes::MS_TOTALMASS] += EMS[e].M;
621}
622}
623
624// Output special points and grids
625if (isSpecialPoint[index_SnowStation]) gatherSpecialPoints(meteoPixel, snowPixel, surfaceFlux);
626fillGrids(ix, iy, meteoPixel, snowPixel, surfaceFlux);
627
628//Restore original water transport scheme that has been changed for ice & firn
629if (land==13 || land==14) {
630sn_cfg.addKey("WATERTRANSPORTMODEL_SNOW", "SnowpackAdvanced", bcu_watertransportmodel_snow);
631sn_cfg.addKey("WATERTRANSPORTMODEL_SOIL", "SnowpackAdvanced", bcu_watertransportmodel_soil);
632}
633//Restore original keys that were modified for glacier pixels
634if (land==14) {
635sn_cfg.addKey("REDUCE_N_ELEMENTS", "SnowpackAdvanced", bcu_reduce_n_elements);
636sn_cfg.addKey("ADJUST_HEIGHT_OF_METEO_VALUES", "SnowpackAdvanced", bcu_adjust_height_of_meteo);
637sn_cfg.addKey("ADJUST_HEIGHT_OF_WIND_VALUE", "SnowpackAdvanced", bcu_adjust_height_of_wind);
638}
639}
640}
641
642void SnowpackInterfaceWorker::grooming(const mio::Date &current_date, const mio::Grid2DObject &grooming_map)
643{
644for(size_t ii = 0; ii < SnowStationsCoord.size(); ++ii) {
645const size_t ix = SnowStationsCoord.at(ii).first;
646const size_t iy = SnowStationsCoord.at(ii).second;
647if (SnowpackInterfaceWorker::skipThisCell(landuse(ix,iy), dem(ix,iy))) continue; //skip nodata cells as well as water bodies, etc
648if (SnowStations[ii]==NULL) continue; //for safety: skipped cells were initialized with NULL
649
650if (grooming_map(ix, iy)==IOUtils::nodata || grooming_map(ix, iy)==0) continue;
651
652if (sn_techsnow.prepare(current_date)) sn_techsnow.preparation(*SnowStations[ii]);
653}
654}
655
656void SnowpackInterfaceWorker::gatherSpecialPoints(const CurrentMeteo& meteoPixel, const SnowStation& snowPixel, const SurfaceFluxes& surfaceFlux)
657{
658meteo_pixel.push_back( meteoPixel );
659snow_pixel.push_back( snowPixel );
660surface_flux.push_back( surfaceFlux );
661}
662
663/**
664 * @brief optimised way to round landuse
665 * @param landuse_dbl is the landuse to round
666 */
667int SnowpackInterfaceWorker::round_landuse(const double& landuse_dbl)
668{
669return (int)(landuse_dbl + 0.0001);
670}
671
672/**
673 * @brief check if a cell should be simulated or skipped
674 * @param landuse_val land use parameter for this pixel
675 * @param dem_val dem altitude for this pixel
676 */
677bool SnowpackInterfaceWorker::skipThisCell(const double& landuse_val, const double& dem_val)
678{
679//determines from the landuse and dem if this cell has to be included in the computation
680//landuse codes are PREVAH codes
681if (landuse_val==IOUtils::nodata) return true;
682
683const int land = (SnowpackInterfaceWorker::round_landuse(landuse_val) - 10000) / 100;
684if (land==9 || land==10 || land==12 || land==16 || land==17 || land>=30) return true; //undefined
685if (land==1) return true;//water
686
687if (dem_val==IOUtils::nodata) return true; //no DEM data
688
689return false;
690}
691
692bool SnowpackInterfaceWorker::is_special(const std::vector< std::pair<size_t,size_t> >& pts_in, const size_t& ix, const size_t& iy)
693{
694for (size_t ii=0; ii<pts_in.size(); ii++) {
695if ((pts_in[ii].first == ix) && (pts_in[ii].second == iy)) return true;
696}
697
698return false;
699}
700
701/**
702 * @brief method called by SnowpackInterface to retrieve all snow_pixels, to read the lateral flow variable
703 * @param ptr_snow_pixel to vector of SnowStations, where the variable SlopeParFlux contains the lateral flow
704 */
705void SnowpackInterfaceWorker::getLateralFlow(std::vector<SnowStation*>& ptr_snow_pixel)
706{
707for(size_t ii = 0; ii < SnowStationsCoord.size(); ++ii) {
708#pragma omp critical (snow_station_lock)
709ptr_snow_pixel.push_back( SnowStations[ii] );
710}
711}
712
713/**
714 * @brief method called by SnowpackInterface to send back the source/sink term for treating lateral flow
715 * @param ptr_snow_pixel to vector of SnowStations, in which the lwc_source variable is set to contain the lateral flow.
716 */
717void SnowpackInterfaceWorker::setLateralFlow(const std::vector<SnowStation*>& ptr_snow_pixel)
718{
719for(size_t ii = 0; ii < SnowStationsCoord.size(); ++ii) {
720if (SnowStations[ii]!=NULL) {
721for(size_t n=0; n<SnowStations[ii]->getNumberOfElements(); n++) {
722SnowStations[ii]->Edata[n].lwc_source = ptr_snow_pixel[ii]->Edata[n].lwc_source;
723}
724}
725}
726}

Archive Download this file

Revision: HEAD