Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/SnowDrift.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/*------------------------------------------------------------------------------------------+
19 | This module is contains the driving routines for a three dimensional numerical model of |
20 | snow drift. It was started in January 00 when Michael had a broken arm and was thinking |
21 | about Betty's surprise birthday party with Fondue which was coming up |
22 +------------------------------------------------------------------------------------------*/
23/********************************************************************************************/
24
25#include <assert.h>
26#include <vector>
27
28#include <alpine3d/snowdrift/SnowDrift.h>
29#include <alpine3d/AlpineMain.h>
30
31#pragma GCC diagnostic push
32#pragma GCC diagnostic ignored "-Weffc++"
33
34using namespace mio;
35using namespace std;
36
37/************************************************************************************/
38
39const double SnowDriftA3D::kinematicViscosityAir = 1.3E-5;//Returns the kinematic viscosity of air (m2 s-1), actually a function of air temperature, The value here is taken from Liston and Sturm (1998)
40const double SnowDriftA3D::USTAR=0.87; //ustar for boundary layer
41const double SnowDriftA3D::molecularWeightofWater = 18.015E-3; //[kg mol-1]
42const double SnowDriftA3D::thermalConductivityofAtm = 0.024;//[J m-1 s-1 K-1] Value taken from Liston and Sturm (1998)
43const double SnowDriftA3D::c_red = 1.0; //Defines red. of saltation conc. in suspension sol.
44const double SnowDriftA3D::grain_size = 0.000680;
45const double SnowDriftA3D::tau_thresh = 0.094; //Original Value by Judith: 0.094
46const double SnowDriftA3D::z0 = 0.01; //Wind Field Z0 - includes larger surface features
47const bool SnowDriftA3D::thresh_snow = true;//Flag to determine whether ustar_thresh is calculated from the Snowpack properties
48
49
50SnowDriftA3D::SnowDriftA3D(const DEMObject& dem, const mio::Config& cfg)
51 : saltation_obj(cfg), auxLayerHeight(0.02), io(cfg), snowpack(NULL), eb(NULL),
52 cH(dem, IOUtils::nodata), sp(dem, IOUtils::nodata), rg(dem, IOUtils::nodata), N3(dem, IOUtils::nodata), rb(dem, IOUtils::nodata),
53 nx(0), ny(0), nz(0), vw(dem, IOUtils::nodata), rh(dem, IOUtils::nodata), ta(dem, IOUtils::nodata), p(dem, IOUtils::nodata),
54 psum(dem, IOUtils::nodata), psum_ph(dem, IOUtils::nodata), STATIONARY(true)
55{
56const string wind_field_string = cfg.get("WINDFIELDS", "Input");
57
58vector<string> TA_interpol;
59cfg.getValues("TA::algorithms", "Interpolations2D", TA_interpol);
60if (TA_interpol.empty())
61throw InvalidArgumentException("No 2D interpolation algorithm defined for TA", AT);
62for (size_t ii=0; ii<TA_interpol.size(); ii++) {
63const string current = TA_interpol[ii];
64if (current!="CST" && current!="AVG" && current!="AVG_LAPSE") {
65if (MPIControl::instance().master())
66cout << "[W] for Snowdrift with sublimation, it is recommended to use CST or AVG or AVG_LAPSE as 2D interpolations for TA\n";
67break;
68}
69}
70
71buildWindFieldsTable(wind_field_string);
72Initialize();
73InitializeFEData();
74}
75
76std::string SnowDriftA3D::getGridsRequirements() const
77{
78return "HS SP RG N3 RB STORE SWE"; //see setSnowSurfaceData() and saltation
79}
80
81void SnowDriftA3D::buildWindFieldsTable(const std::string& wind_field_string)
82{//parses the given string containing the wind fields and number of steps, for example: "NW 3 SE 5 N 6"
83//the result is written into a vector of wind_field structures, wind_fields
84int previous_time_step = 0, steps, fields_count=0;
85std::istringstream iss(wind_field_string, std::istringstream::in);
86struct WIND_FIELD curr_field;
87
88std::string word;
89while ( iss >> word ) {
90curr_field.wind = word;
91if (!(iss>>word)) {
92throw InvalidArgumentException("Missing number of time steps for wind field "+curr_field.wind, AT);
93}
94IOUtils::convertString(steps, word);
95curr_field.start_step = previous_time_step;
96previous_time_step += steps;
97fields_count ++;
98wind_fields.push_back(curr_field);
99}
100wind_field_index = -1;
101cout << "[I] Snowdrift will use " << fields_count << " wind situations" << endl;
102}
103
104bool SnowDriftA3D::isNewWindField(const unsigned int current_step)
105{//this function returns true if a new wind field should be read
106const unsigned int next_index = wind_field_index+1;
107
108if (wind_field_index==-1) return true;
109
110if ((next_index<wind_fields.size()) && (current_step >= wind_fields[next_index].start_step)) {
111return true;
112} else {
113return false;
114}
115}
116
117/**
118 * @brief Initialize for snowdrift
119 * Allocate saltation variable, snow height and new snow mass per bottom element
120 */
121void SnowDriftA3D::Initialize()
122{
123Grid3DObject z_readMatr;
124io.read3DGrid(z_readMatr, wind_fields[0].wind+":DEM");
125z_readMatr.llcorner = ta.llcorner; //most probably the user did not properly specify the ARPS_XCOORDS and ARPS_YCOORDS
126
127nx = z_readMatr.getNx();
128ny = z_readMatr.getNy();
129nz = z_readMatr.getNz();
130
131/* Allocate some memory and Initialize */
132/* Before entering the time integration loop allocate the saltation variable,
133the snow height and new snow mass per bottom element */
134saltation.resize(nx,ny, 0.);
135c_salt.resize(nx,ny);
136
137mns.set(nx, ny, ta.cellsize, ta.llcorner, 0.);
138mns_subl.resize(nx,ny);
139mns_nosubl.resize(nx,ny);
140dif_mns_subl.resize(nx,ny);
141flux_x.resize( nx*ny );
142flux_y.resize( nx*ny );
143flux_z.resize( nx*ny );
144flux_x_subl.resize( nx*ny );
145flux_y_subl.resize( nx*ny );
146flux_z_subl.resize( nx*ny );
147
148InitializeNodes(z_readMatr);
149ConstructElements();
150}
151
152void SnowDriftA3D::Destroy() {}
153
154double SnowDriftA3D::getTiming() const
155{
156return timer.getElapsed();
157}
158
159SnowDriftA3D::~SnowDriftA3D()
160{
161Destroy();
162}
163
164/**CompleteNodes
165 * @brief Initialize concentration according to the rain rate (uniform)
166 * Now calculate the missing parameters
167 * Slope in direction of the wind
168 */
169void SnowDriftA3D::CompleteNodes()
170{
171if ( new_wind_status == true ) {
172
173/* The turbulent settling velocity for the whole domain is now set to WS0 (formerly it
174contained the */
175/* Look at subscale turbulence effect */
176nodes_wstar.grid3D.resize(nx,ny,nz, WS0);
177
178// LH_CHANGE BEGIN: re-definition of the arps-mesh; Adapted by ML on August 12 2006
179cout << "[i] Snowdrift: updating the mesh..."<<endl;
180
181// the second arps layer (topography) becomes the bottom (ie boundary)
182// layer of the suspension grid, calculate the slope first
183for (size_t ii=0; ii<nx; ii++){
184for (size_t jj=0; jj<ny; jj++){
185nodes_slope(ii,jj,0) = (nodes_v(ii,jj,1)*atan(-nodes_sy(ii,jj,1)) + nodes_u(ii,jj,1)*atan(-nodes_sx(ii,jj,1)))
186/sqrt(Optim::pow2(nodes_v(ii,jj,1)) + Optim::pow2(nodes_u(ii,jj,1)));
187nodes_slope(ii,jj,1) = nodes_slope(ii,jj,0);
188nodes_u(ii,jj,0) = 0.;
189nodes_v(ii,jj,0) = 0.;
190nodes_w(ii,jj,0) = 0.;
191nodes_wstar(ii,jj,0) = nodes_wstar(ii,jj,1);
192nodes_e(ii,jj,0) = nodes_e(ii,jj,1);
193}
194}
195
196//an additional suspension grid layer of nodes is added between the
197//first and second arps layer
198
199//the new layer is located between the first and second
200//arps layer, adjusted by the factor auxLayerHeight. The x,y
201//coordinates remain the same as long as the mesh remains
202//regular in these directions
203
204for (size_t ii=0; ii<nx; ii++){
205for (size_t jj=0; jj<ny; jj++){
206const double salt_height = nodes_z(ii,jj,1) - nodes_z(ii,jj,0);
207
208//the direction of the wind field of a node in the new layer is
209//assumed to be equal to the node lying in the next higher layer above.
210//and the magnitude is scaled by a factor according to a logarithmic wind profile
211const double fac = log( salt_height/z0 ) / log( salt_height/(auxLayerHeight * z0) );
212
213nodes_u(ii,jj,1) = fac*nodes_u(ii,jj,2);
214nodes_v(ii,jj,1) = fac*nodes_v(ii,jj,2);
215nodes_w(ii,jj,1) = fac*nodes_w(ii,jj,2);
216
217//the turbulent kinetic energy of the new layer is given by a
218//weighted average
219nodes_e(ii,jj,1) = (1.-auxLayerHeight) *nodes_e(ii,jj,0) + auxLayerHeight * nodes_e(ii,jj,2);
220//nodes_slope(ii,jj,1) = nodes_slope(ii,jj,2);
221}
222}
223} // End if new wind
224}
225
226/**
227 * @brief Initialize nodes
228 * Each (3D) node receives its (x,y,z), slope, aspect, normals
229 */
230void SnowDriftA3D::InitializeNodes(const mio::Grid3DObject& z_readMatr)
231{
232//we create a DEMObject using the Corripio slope algorithm
233const double cellsize = z_readMatr.cellsize;
234
235nodes_z = z_readMatr; //full copy
236nodes_x.set(z_readMatr, 0.);
237nodes_y.set(z_readMatr, 0.);
238nodes_sx.set(z_readMatr, 0.);
239nodes_sy.set(z_readMatr, 0.);
240nodes_e.set(z_readMatr, 0.);
241nodes_c.set(z_readMatr, 0.);
242nodes_Tair.set(z_readMatr, 0.);
243nodes_slope.set(z_readMatr, 0.);
244
245if (SUBLIMATION) {
246nodes_RH.set(z_readMatr, 0.); //relative humidity
247nodes_q.set(z_readMatr, 0.);
248nodes_Subl.set(z_readMatr, 0.);
249
250nodes_tmp_c.set(z_readMatr, 0.);
251nodes_Tair_ini.set(z_readMatr, 0.);
252nodes_q_ini.set(z_readMatr, 0.);
253nodes_Subl_ini.set(z_readMatr, 0.);
254}
255
256for (unsigned int kk=0; kk<nz; kk++){
257for (unsigned int jj=0;jj<ny;jj++){
258for (unsigned int ii=0;ii<nx;ii++){
259nodes_x(ii,jj,kk) = ii*cellsize;
260nodes_y(ii,jj,kk) = jj*cellsize;
261}
262}
263}
264
265DEMObject curr_layer(ta, false, DEMObject::CORR); //take the geolocalization from ta
266curr_layer.setUpdatePpt((DEMObject::update_type)(DEMObject::SLOPE | DEMObject::NORMAL));
267for (unsigned int kk=0; kk<4; kk++){
268nodes_z.extractLayer(kk, curr_layer);
269curr_layer.update(); //compute the slopes, normals, min/max, etc only for the bottom 3 layers
270for (unsigned int jj=0;jj<ny;jj++){
271for (unsigned int ii=0;ii<nx;ii++){
272nodes_sx(ii,jj,kk) = curr_layer.Nx(ii,jj);
273nodes_sy(ii,jj,kk) = curr_layer.Ny(ii,jj);
274}
275}
276}
277
278// LH_CHANGE BEGIN: re-definition of the arps-mesh
279cout <<"[i] Snowdrift adding artificial layer..."<<endl;
280
281// the second arps layer (topography) becomes the bottom (ie boundary)
282// layer of the suspension grid
283for (unsigned int ii=0;ii<nx;ii++){
284for (unsigned int jj=0;jj<ny;jj++){
285nodes_z(ii,jj,0) = z_readMatr(ii,jj,1);
286nodes_x(ii,jj,0) = nodes_x(ii,jj,1);
287nodes_y(ii,jj,0) = nodes_y(ii,jj,1);
288nodes_sx(ii,jj,0) = nodes_sx(ii,jj,1);
289nodes_sy(ii,jj,0) = nodes_sy(ii,jj,1);
290}
291}
292//an additional suspension grid layer of nodes is added between the
293//first and second arps layer
294//the new layer is located between the first and second
295//arps layer, adjusted by the factor auxLayerHeight. The x,y
296//coordinates remain the same as long as the mesh remains
297//regular in these directions
298
299for (unsigned int ii=0;ii<nx;ii++){
300 for (unsigned int jj=0;jj<ny;jj++){
301const double salt_height = auxLayerHeight*(nodes_z(ii,jj,2) - nodes_z(ii,jj,0));
302if (salt_height<0.) std::cout << "[E] Invalid height ARPS data at (" << ii << "," << jj << ") = " << salt_height << "\n";
303nodes_z(ii,jj,1) = nodes_z(ii,jj,0)+salt_height;
304nodes_x(ii,jj,1) = nodes_x(ii,jj,2);
305nodes_y(ii,jj,1) = nodes_y(ii,jj,2);
306nodes_sx(ii,jj,1) = nodes_sx(ii,jj,2);
307nodes_sy(ii,jj,1) = nodes_sy(ii,jj,2);
308 }
309}
310}
311
312void SnowDriftA3D::ConstructElements()
313{
314//arrays for mapping local node indices onto global ones
315dofMap.resize( (nx-1)*(ny-1)*(nz-1) + 1, 8);
316/*
317* the dofMap matrix is used together with c and Psi, i.e. with vectors
318* that contain only the degrees of freedom; for building up dofMap we take
319* all elements of the domain but the application does only make
320* sense for nodes that are on the interior of the domain; we go from
321* 1 to (nx-1)*(ny-1)*(nz-1) elements and assign them their nodes
322* (from 1 to (nx-2)*(ny-2)*(nz-2) nodes)
323*/
324for (unsigned int kk = 0; kk < (nz-1); kk++ ) {
325for (unsigned int jj = 0; jj < (ny-1); jj++ ) {
326for (unsigned int ii = 1; ii <= (nx-1); ii++ ) {
327const int elem = kk*(nx-1)*(ny-1) + jj*(nx-1) + ii;
328
329dofMap(elem,0) = (jj-1)*(nx-2) + ii-1 + (kk-1)*(nx-2)*(ny-2);
330dofMap(elem,1) = (jj-1)*(nx-2) + (ii+0) + (kk-1)*(nx-2)*(ny-2);
331dofMap(elem,2) = (jj)*(nx-2) + (ii+0) + (kk-1)*(nx-2)*(ny-2);
332dofMap(elem,3) = (jj)*(nx-2) + ii-1 + (kk-1)*(nx-2)*(ny-2);
333dofMap(elem,4) = (jj-1)*(nx-2) + ii-1 + kk*(nx-2)*(ny-2);
334dofMap(elem,5) = (jj-1)*(nx-2) + (ii+0) + kk*(nx-2)*(ny-2);
335dofMap(elem,6) = (jj)*(nx-2) + (ii+0) + kk*(nx-2)*(ny-2);
336dofMap(elem,7) = (jj)*(nx-2) + ii-1 + kk*(nx-2)*(ny-2);
337}
338}
339}
340
341
342/*
343* nodeMap
344* this is the nodeMap matrix with all the elements from 1 to
345* (nx-1)*(ny-1)*(nz-1) and a mapping on all nodes of the domain note
346* that this approach is good for the enumeration of the real ARPS
347* mesh with the underground layer and is also good shit for the mesh
348* used for the calculation, i.e. without underground but with the
349* artificial nodes layer
350*/
351nodeMap.resize( (nx-1)*(ny-1)*(nz-1), 8);
352int element = 0;
353for (unsigned int iz=0; iz<(nz-1); iz++) {
354for (unsigned int iy=0; iy<(ny-1); iy++) {
355for (unsigned int ix=0; ix<(nx-1); ix++) {
356nodeMap(element,0) = ix + iy*nx + iz*nx*ny;
357nodeMap(element,1) = nodeMap(element,0) + 1;
358nodeMap(element,2) = nodeMap(element,0) + nx+1;
359nodeMap(element,3) = nodeMap(element,0) + nx;
360nodeMap(element,4) = nodeMap(element,0) + nx*ny;
361nodeMap(element,5) = nodeMap(element,4) + 1;
362nodeMap(element,6) = nodeMap(element,4) + nx+1;
363nodeMap(element,7) = nodeMap(element,4) + nx;
364element++;
365}
366}
367}
368
369
370elems.resize((nx-1)*(ny-1)*(nz-1), 8);
371int nelems=0;
372for (unsigned int iz=0; iz<(nz-1); iz++) {
373for (unsigned int iy=0; iy<(ny-1); iy++) {
374for (unsigned int ix=0; ix<(nx-1); ix++) {
375elems(nelems,0) = ix + iy*nx + iz*nx*ny;
376elems(nelems,4) = elems[nelems][0] + nx*ny;
377elems(nelems,1) = elems[nelems][0] + 1;
378elems(nelems,3) = elems[nelems][0] + nx;
379elems(nelems,2) = elems[nelems][0] + nx+1;
380elems(nelems,5) = elems[nelems][4] + 1;
381elems(nelems,7) = elems[nelems][4] + nx;
382elems(nelems,6) = elems[nelems][4] + nx+1;
383++nelems;
384}
385}
386}
387}
388
389
390void SnowDriftA3D::SnowMassChange(bool setbound, const mio::Date& calcDate)
391{
392//LH_DEBUG BEGIN
393static int timestepp=0;
394timestepp++;
395
396//loop over all interior nodes, boundary nodes are set to zero
397//subsequently
398for (unsigned int iy = 1; iy < ny-1; iy++) {
399for (unsigned int ix = 1; ix < nx-1; ix++) {
400dif_mns_subl(ix,iy)=0.; //reset
401if (!SALTATION) {
402mns(ix,iy) = 0.;
403continue; //saltation is not computed, only initiliazing
404}
405
406// second ARPS layer is topography
407// Now first layer is on the ground and has a ground parallel velocity as Henning has tested
408// n = iy*nx + ix;
409
410// Enumerate the elements and nodes used to calculate saltation;
411// units of saltation[i][j] are kg/m/s
412double salt[9];
413salt[0] = salt[4] = saltation(ix-1,iy);
414salt[1] = saltation(ix,iy-1);
415salt[2] = saltation(ix+1,iy);
416salt[3] = saltation(ix,iy+1);
417salt[5] = saltation(ix-1,iy-1);
418salt[6] = saltation(ix+1,iy-1);
419salt[7] = saltation(ix+1,iy+1);
420salt[8] = saltation(ix-1,iy+1);
421
422Array2D<int> grid_ixiy(10,2);//tested Array2D<int> grid_ixiy(2,10);
423
424grid_ixiy(0,0)=ix-1;
425grid_ixiy(0,1)=iy;
426
427grid_ixiy(1,0)=ix;
428grid_ixiy(1,1)=iy-1;
429
430grid_ixiy(2,0)=ix+1;
431grid_ixiy(2,1)=iy;
432
433grid_ixiy(3,0)=ix;
434grid_ixiy(3,1)=iy+1;
435
436grid_ixiy(4,0)=ix-1;
437grid_ixiy(4,1)=iy;
438
439grid_ixiy(5,0)=ix;
440grid_ixiy(5,1)=iy-1;
441
442grid_ixiy(6,0)=ix-1;
443grid_ixiy(6,1)=iy-1;
444
445grid_ixiy(7,0)=ix+1;
446grid_ixiy(7,1)=iy-1;
447
448grid_ixiy(8,0)=ix+1;
449grid_ixiy(8,1)=iy+1;
450
451grid_ixiy(9,0)=ix-1;
452grid_ixiy(9,1)=iy+1;
453
454const double DX=nodes_x.grid3D(grid_ixiy(2,0),grid_ixiy(2,1),1)-nodes_x.grid3D(grid_ixiy(1,0),grid_ixiy(1,1),1);
455const double DY=nodes_y.grid3D(grid_ixiy(2,0),grid_ixiy(2,1),1)-nodes_y.grid3D(grid_ixiy(1,0),grid_ixiy(1,1),1);
456const double DZ=nodes_z.grid3D(grid_ixiy(2,0),grid_ixiy(2,1),1)-nodes_z.grid3D(grid_ixiy(1,0),grid_ixiy(1,1),1);
457
458const double DX1=nodes_x.grid3D(grid_ixiy(0,0),grid_ixiy(0,1),1)-nodes_x.grid3D(grid_ixiy(1,0),grid_ixiy(1,1),1);
459const double DY1=nodes_y.grid3D(grid_ixiy(0,0),grid_ixiy(0,1),1)-nodes_y.grid3D(grid_ixiy(1,0),grid_ixiy(1,1),1);
460const double DZ1=nodes_z.grid3D(grid_ixiy(0,0),grid_ixiy(0,1),1)-nodes_z.grid3D(grid_ixiy(1,0),grid_ixiy(1,1),1);
461
462double f_cross[3];
463f_cross[0] = DY*DZ1-DZ*DY1; /* cross product of DX(vector) and DX1(vector) */
464f_cross[1] = DZ*DX1-DX*DZ1;
465f_cross[2] = DX*DY1-DY*DX1;
466
467const double area = sqrt( Optim::pow2(f_cross[0]) + Optim::pow2(f_cross[1]) + Optim::pow2(f_cross[2]) );
468
469/* area of the element is the norm of the cross product */
470double salt_flux = 0.0; /* The four contributions to erosion / deposition at a node */
471
472//now the edge loop in order to calculate the contribution
473//of each edge to the final salt_flux of the element
474double salt_flux_frac[4];
475
476double fac1[4];
477for (size_t i = 0; i < 4; i++) {
478//4ptdiv
479const double u = 0.25 * ( nodes_u.grid3D(ix,iy,1) + nodes_u.grid3D(grid_ixiy(i,0),grid_ixiy(i,1),1) + nodes_u.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) + nodes_u.grid3D(grid_ixiy(i+6,0),grid_ixiy(i+6,1),1));
480const double v = 0.25 * ( nodes_v.grid3D(ix,iy,1) + nodes_v.grid3D(grid_ixiy(i,0),grid_ixiy(i,1),1) + nodes_v.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) + nodes_v.grid3D(grid_ixiy(i+6,0),grid_ixiy(i+6,1),1));
481const double w = 0.25 * ( nodes_w.grid3D(ix,iy,1) + nodes_w.grid3D(grid_ixiy(i,0),grid_ixiy(i,1),1) + nodes_w.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) + nodes_w.grid3D(grid_ixiy(i+6,0),grid_ixiy(i+6,1),1));
482
483// vector from node i to i+1
484const double dx = ( nodes_x.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) - nodes_x.grid3D(grid_ixiy(i,0),grid_ixiy(i,1),1));
485const double dy = ( nodes_y.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) - nodes_y.grid3D(grid_ixiy(i,0),grid_ixiy(i,1),1));
486const double dz = ( nodes_z.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) - nodes_z.grid3D(grid_ixiy(i,0),grid_ixiy(i,1),1));
487
488// only valid if grid is rectangular, take vector from
489// node i+1 to node i+2 as the normal vector on (dx,dy,dz)
490const double dx1 = ( nodes_x.grid3D(grid_ixiy(i+2,0),grid_ixiy(i+2,1),1) - nodes_x.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) );
491const double dy1 = ( nodes_y.grid3D(grid_ixiy(i+2,0),grid_ixiy(i+2,1),1) - nodes_y.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) );
492const double dz1 = ( nodes_z.grid3D(grid_ixiy(i+2,0),grid_ixiy(i+2,1),1) - nodes_z.grid3D(grid_ixiy(i+1,0),grid_ixiy(i+1,1),1) );
493
494// if there is no wind, set fac1 to zero
495// LH_REMARK: testing floats on equality? dangerous!
496if ( (u == 0.0) && (v == 0.0) && (w==0.0) ) {
497printf("Scheiss\n");
498 fac1[i] = 0;
499} else {
500//fac1=velocity(unitvec.) * lineunitvector that is perp.to
501//dx(vec)*abs(dx(vector))/area of the element;
502//the units of fac1 are 1/m, or equivalently: *m/m2 makes more physical sense
503fac1[i] = (u * dx1 + v * dy1 + w * dz1)
504/ sqrt( Optim::pow2(u) + Optim::pow2(v) + Optim::pow2(w) )
505/ sqrt( Optim::pow2(dx1) + Optim::pow2(dy1) + Optim::pow2(dz1) )
506/ area * sqrt( Optim::pow2(dx) + Optim::pow2(dy) + Optim::pow2(dz) );
507}
508salt_flux_frac[i] = fac1[i] * 1/4* (salt[i]+salt[i+1]+salt[i+5]+saltation(ix,iy)); /*flux over 1 of the four edges*/
509
510if (salt_flux_frac[i] > 0 && (salt[i+5] == 0)){
511// point at other side doesn't have saltation, block incoming flux
512salt_flux_frac[i] = 0;
513}
514if (salt_flux_frac[i] < 0 && (saltation(ix,iy) == 0)){
515// erosion necessary but there was no snow, block outgoing flux
516salt_flux_frac[i] = 0;
517}
518
519salt_flux+=salt_flux_frac[i];/*sum over the four edges*/
520
521// here we calculate the saltation flux on edge i by multiplication
522// of the saltation by the lenght of the effective edge and by division
523// of the area of the element, therefore we obtain a flux in units of kg/m2/s
524
525} // end edge loop
526
527
528// scalar product of flux and normal unit vector to the surface, -ve
529// sign because the unit vector is in +ve z direction; units of flux_i[j]
530// is kg/m2/s , therefore also for diff_flux because we multiplicate
531// by unit vector f that is dimensionless ;
532/* The four contributions to erosion / deposition at a node */
533double diff_flux = -( flux_x[nx*iy + ix] * f_cross[0]
534+ flux_y[nx*iy + ix] * f_cross[1]
535+ flux_z[nx*iy + ix] * f_cross[2] ) / area;
536
537const double diff_flux_subl = (SUBLIMATION)? -( flux_x_subl[nx*iy + ix] * f_cross[0]
538+ flux_y_subl[nx*iy + ix] * f_cross[1] + flux_z_subl[nx*iy + ix] * f_cross[2] ) / area : 0.;
539
540mns(ix,iy) = dt_main * ( salt_flux + diff_flux );
541if (SUBLIMATION){
542mns_subl(ix,iy) = dt_main * ( salt_flux + diff_flux_subl);
543//cout<<"mns (nosubl)= " <<mns(ix,iy)<<" mns (subl)=" << mns_subl(ix,iy)<<endl;
544}
545/*the total flux is the sum of saltation and diff_flux*/
546
547//this test to limit crazy drift is normally done in SnowInterface, do this here to make sure to be able to compare sims with and without subl
548if (mns(ix,iy) > 37.) {
549//printf("Crazy drift: ix %d, iy %d mm/h %f\n", ix, iy, mns(ix,iy));
550mns(ix,iy) = 37.;
551}
552if (mns(ix,iy) < -37.) {
553//printf("Crazy drift: ix %d, iy %d mm/h %f\n", ix, iy, mns(ix,iy));
554mns(ix,iy) = -37.;
555}
556
557if (SUBLIMATION){
558if (mns_subl(ix,iy) > 37.) {
559printf("Crazy drift (subl): ix %d, iy %d mm/h %f\n", ix, iy, mns_subl(ix,iy));
560mns_subl(ix,iy) = 37.;
561}
562if (mns_subl(ix,iy) < -37.) {
563printf("Crazy drift (subl): ix %d, iy %d mm/h %f\n", ix, iy, mns_subl(ix,iy));
564mns_subl(ix,iy) = -37.;
565}
566
567dif_mns_subl(ix,iy)=mns(ix,iy)-mns_subl(ix,iy);
568}else{
569dif_mns_subl(ix,iy)=0.;
570}
571
572if (SUBLIMATION_OUTPUT){
573mns_nosubl(ix,iy)=mns(ix,iy);//just to save this in file
574}
575
576if (SUBLIMATION){
577//take value with sublimation as standard
578mns(ix,iy)=mns_subl(ix,iy);
579diff_flux=diff_flux_subl;
580}
581
582 //if ( (diff_flux < 0.0) && (saltation(ix,iy) == 0.0) && (salt_flux == 0.0)) {
583// mns(ix,iy) = std::max(0.0,mns(ix,iy));
584// }
585
586#if WRITE_DRIFT_FLUXES
587const string fname = "../output/" + calcDate.toString(Date::NUM) + ".dep";
588FILE *difFile = fopen(fname.c_str(), "a");
589if (difFile == NULL)
590throw AccessException("Can not open file '"+fname+"'", AT);
591fprintf(difFile, "%s %d %d %f %f %f\n", calcDate.toString(Date::ISO).c_str(), ix, iy ,dt_main*salt_flux, dt_main*diff_flux, mns(ix,iy));
592fclose( difFile );
593#endif
594if (fabs(mns(ix,iy)) > 1000.) {
595cout<<"Halt Salt:"<<ix<<" "<<iy<<" "<<salt_flux<<" "<<diff_flux<<endl;
596}
597}
598}/* End node loop */
599
600
601if (SUBLIMATION_OUTPUT){
602//Allows to obtain sublimation amount in 1 simulation > used for season simulation. NOTE: needs to be corrected for slope when comparing to SWE-grid
603const string fname= string("../results/")+ calcDate.toString(Date::NUM)+ string(".dif_subl");
604FILE *difFile = fopen(fname.c_str(), "w");
605if (difFile == NULL) {
606cout <<"Cannot open file "<<fname.c_str()<<endl;
607return;
608}
609for (unsigned int iy = 0; iy < ny; iy++){
610for (unsigned int ix = 0;ix < nx; ix++){
611fprintf(difFile, " %f", dif_mns_subl(ix,iy));
612}
613fprintf(difFile,"\n");
614}
615fclose(difFile);
616
617string fname2= string("../results/")+ calcDate.toString(Date::NUM)+ string(".mns_nosubl");
618FILE *mnsFile;
619mnsFile=fopen(fname2.c_str(), "w");
620if (mnsFile == NULL) {
621cout <<"Cannot open file "<<fname2.c_str()<<endl;
622return;
623}
624
625for (unsigned int iy = 0; iy < ny; iy++){
626for (unsigned int ix = 0;ix < nx; ix++){
627fprintf(mnsFile, " %f", mns_nosubl(ix,iy));
628}
629fprintf(mnsFile,"\n");
630}
631fclose(mnsFile);
632}
633
634// Now set border values to zero
635if ( setbound ) {
636for (unsigned int iy = 0 ; iy < ny; iy++ ) {
637mns(0,iy) = 0.0;
638mns(nx-1,iy) = 0.0;
639mns(1,iy) = 0.0;
640mns(nx-2,iy) = 0.0;
641}
642for (unsigned int ix = 0; ix < nx; ix++ ) {
643mns(ix,0) = 0.0;
644mns(ix,ny-1) = 0.0;
645mns(ix,1) = 0.0;
646mns(ix,ny-2) = 0.0;
647}
648}
649}
650
651/**
652 * @brief Main: Calls the essential routines
653 * @param calcDate date of current time step
654 */
655void SnowDriftA3D::Compute(const Date& calcDate)
656{
657timer.restart();
658
659const double max_hs = cH.grid2D.getMax();
660const double max_psum = psum.grid2D.getMax();
661const double min_psum_ph = psum_ph.grid2D.getMin();
662if (max_hs==0. && (max_psum==0. || min_psum_ph>=1.)) {
663cout<<"[i] SnowDrift: no snow on the ground & no solid precipitation\n";
664mns=0.;
665if (snowpack!=NULL) {
666snowpack->setSnowMassChange(mns, calcDate);
667}
668return;
669}
670
671/* Calculate the Saltation Fluxes */
672if (SALTATION) {
673cout<<"[i] SnowDrift starting Saltation\n";
674compSaltation(true);
675}
676
677/* Calculate the suspension drift for the wind field until stationary */
678cout<<"[i] SnowDrift starting Suspension\n";
679// DEBUG("***Checksum BEFORE Suspension: nodes=%lf, nodes[c]=%lf", checksum(nodes), checksum_c(nodes_c));
680Suspension();
681// DEBUG("***Checksum after Suspension: nodes=%lf, nodes[c]=%lf", checksum(nodes), checksum_c(nodes_c));
682
683if (FIELD3D_OUTPUT){
684 debugOutputs(calcDate, string("../results/"), OUT_CONC);
685 }
686
687SnowMassChange(true, calcDate);
688
689std::cout << "[i] SnowDrift simulation done for " << calcDate.toString(Date::ISO) << "\n";
690
691if (snowpack!=NULL) {
692snowpack->setSnowMassChange(mns, calcDate);
693}
694
695//DEBUG("Checksum SnowDrift: nodes=%lf, saltation=%lf, c_salt=%lf", checksum_c(nodes_c), checksum(saltation), checksum(c_salt));
696timer.stop();
697
698} /* End SnowDrift */
699
700void SnowDriftA3D::GetTResults(double outtime_v[15], double outtime_tau[15], double outtime_salt[15], double outtime_diff[15])
701{
702const size_t sz = 15*sizeof(double);
703memcpy(outtime_v, time_v, sz);
704memcpy(outtime_tau, time_tau, sz);
705memcpy(outtime_salt, time_salt, sz);
706memcpy(outtime_diff, time_diff, sz);
707}
708
709void SnowDriftA3D::setSnowPack(SnowpackInterface &mysnowpack)
710{
711snowpack=&mysnowpack;
712}
713
714void SnowDriftA3D::setEnergyBalance(EnergyBalance &myeb)
715{
716eb=&myeb;
717}
718
719void SnowDriftA3D::setSnowSurfaceData(const mio::Grid2DObject& cH_in, const mio::Grid2DObject& sp_in, const mio::Grid2DObject& rg_in,
720 const mio::Grid2DObject& N3_in, const mio::Grid2DObject& rb_in)
721{
722cH = cH_in;
723sp = sp_in;
724rg = rg_in;
725N3 = N3_in;
726rb = rb_in;
727}
728
729void SnowDriftA3D::debugOutputs(const Date& calcDate, const std::string& outpath, const DRIFT_OUTPUT& output_type)
730{
731const std::string ext = (output_type==OUT_CONC)? ".ctn" : ".sub";
732const std::string fname( outpath + calcDate.toString(Date::NUM) + ext );
733
734writeOutput(fname); //write output file of snowdrift
735}
736
737
738//HACK: this should be done by MeteoIO
739/**
740 * @brief Write output
741 * Writes the values of several nodes-fields
742 */
743void SnowDriftA3D::writeOutput(const std::string& fname)
744{
745FILE *logfile = fopen(fname.c_str(), "w");
746if (logfile == NULL) {
747cout<<"Cannot open file"<<fname.c_str()<<endl;
748return;
749}
750fprintf(logfile, "#1:j 2:ix 3:iy 4:iz \t 5:z 6:c");
751if (SUBLIMATION_OUTPUT) fprintf(logfile, " \t 7:sublimation 8:RH 9:q 10:Ta");
752fprintf(logfile, " \t 11:u 12:v 13:w \n");
753
754for (unsigned int iy=0; iy<ny; iy++) {
755for (unsigned int ix=0; ix<nx; ix++) {
756 //compute int dz c
757 double intc=0;
758 for (unsigned int iz=0; iz<nz-1; iz++) {
759const int j = iz*nx*ny + iy*nx + ix;
760 intc += c[j]*(nodes_z.grid3D(ix,iy,iz+1)-nodes_z.grid3D(ix,iy,iz));
761 }
762 intc /= (nodes_z.grid3D(ix,iy,nz-1)-nodes_z.grid3D(ix,iy,0));
763
764 for (unsigned int iz=0; iz<nz; iz++) {
765const int j = iz*nx*ny + iy*nx + ix;
766fprintf(logfile, "%d %d %d %d \t %f %f", j, ix, iy, iz, nodes_z.grid3D(ix,iy,iz), c[j]*1000);
767
768if (SUBLIMATION_OUTPUT)
769fprintf(logfile, " \t %f %f %f %f",
770nodes_Subl.grid3D(ix,iy,iz)*1000, nodes_RH.grid3D(ix,iy,iz), nodes_q.grid3D(ix,iy,iz)*1000, nodes_Tair.grid3D(ix,iy,iz));
771
772fprintf(logfile, " \t %f %f %f\n", nodes_u.grid3D(ix,iy,iz), nodes_v.grid3D(ix,iy,iz), nodes_w.grid3D(ix,iy,iz));
773 }
774}
775}
776fclose (logfile);
777}
778
779/**
780* @brief Sets the required meteo fields
781*/
782void SnowDriftA3D::setMeteo (const unsigned int& steps, const Grid2DObject& new_psum, const mio::Grid2DObject& new_psum_ph, const Grid2DObject& new_p, const Grid2DObject& /*new_vw*/,
783 const Grid2DObject& new_rh, const Grid2DObject& new_ta, const Grid2DObject& new_ilwr, const mio::Date& calcDate,
784 const std::vector<mio::MeteoData>& vecMeteo)
785{
786if (vecMeteo.empty())
787throw mio::NoDataException("No meteo data!", AT);
788
789// find the first (often the only one) MeteoData having iswr, ta and ilwr
790// it should exist, otherwise the checkInputsRequirements method would have failed.
791bool found = false;
792for (size_t i = 0; i < vecMeteo.size() ; ++i) {
793ta_1D = vecMeteo[i](MeteoData::TA);
794station_altitude = vecMeteo[i].meta.position.getAltitude();
795if (ta_1D != IOUtils::nodata && station_altitude!=IOUtils::nodata) {
796found = true;
797break;
798}
799}
800if (!found)
801throw mio::NoDataException("No TA station could be found at "+vecMeteo[0].date.toString(mio::Date::ISO), AT);
802
803//Common meteo fields
804psum = new_psum;
805psum_ph = new_psum_ph;
806rh = new_rh;
807ta = new_ta;
808p = new_p;
809new_wind_status = isNewWindField(steps);
810
811if (new_wind_status) {
812wind_field_index++;
813const std::string filename( wind_fields[wind_field_index].wind );
814io.read3DGrid(nodes_u, filename+":u");
815io.read3DGrid(nodes_v, filename+":v");
816io.read3DGrid(nodes_w, filename+":w");
817if (READK) io.read3DGrid(nodes_K, filename+":kmh");
818cout <<"[i] Snowdrift: ARPS wind field successfully read"<<endl;
819}
820
821mio::Grid2DObject dw(vw, IOUtils::nodata);// dw field with vw as template for dimensions
822//adjust the meteo fields that depend on the 3D wind field
823for (size_t jj=0; jj<ny; jj++) {
824for (size_t ii=0; ii<nx; ii++) {
825vw.grid2D(ii,jj) = Optim::fastSqrt_Q3( Optim::pow2(nodes_u.grid3D(ii,jj,2)) + Optim::pow2(nodes_v.grid3D(ii,jj,2)) ); //Third layer is first layer in the air
826dw.grid2D(ii,jj) = (atan2(nodes_u.grid3D(ii,jj,2), nodes_v.grid3D(ii,jj,2))) * mio::Cst::to_deg;
827}
828}
829
830CompleteNodes();
831if (SUBLIMATION) initializeTRH();
832
833//TODO: feedback mecanism: make it more general!
834if (snowpack!=NULL) snowpack->setMeteo(psum, psum_ph, vw, dw, rh, ta, calcDate);
835if (eb!=NULL) eb->setMeteo(new_ilwr, ta, rh, p, calcDate);
836}
837
838/**
839* @brief Initialize RH, T and q
840* Initialize fields of humidity and temperature, necessary for sublimation calculations
841*/
842void SnowDriftA3D::initializeTRH()
843{
844//start with first level
845for (unsigned int ix=0; ix<nx; ix++){
846for (unsigned int iy=0; iy<ny; iy++){
847nodes_RH.grid3D(ix,iy,0) = rh.grid2D(ix,iy); //constant field, single measurement at Wan3
848nodes_Tair.grid3D(ix,iy,0) = ta_1D;
849nodes_q.grid3D(ix,iy,0) = Atmosphere::relToSpecHumidity(nodes_z.grid3D(ix,iy,0), nodes_Tair.grid3D(ix,iy,0),nodes_RH.grid3D(ix,iy,0));
850nodes_Tair.grid3D(ix,iy,0) = ta_1D-(nodes_z.grid3D(ix,iy,0) - station_altitude)*Cst::dry_adiabatique_lapse_rate;
851nodes_RH.grid3D(ix,iy,0) = RH_from_q(nodes_Tair.grid3D(ix,iy,0),nodes_q.grid3D(ix,iy,0), nodes_z.grid3D(ix,iy,0));
852}
853}
854
855//adjust rest of domain
856for (unsigned int iz=1; iz<nz; iz++){
857for (unsigned int iy=0; iy<ny; iy++){
858for (unsigned int ix=0;ix<nx; ix++){
859nodes_q.grid3D(ix,iy,iz)=nodes_q.grid3D(ix,iy,iz-1);
860nodes_Tair.grid3D(ix,iy,iz)=ta_1D-(nodes_z.grid3D(ix,iy,iz) - station_altitude)*Cst::dry_adiabatique_lapse_rate;
861nodes_RH.grid3D(ix,iy,iz)=RH_from_q(nodes_Tair.grid3D(ix,iy,iz),nodes_q.grid3D(ix,iy,iz), nodes_z.grid3D(ix,iy,iz));
862}
863}
864}
865nodes_Tair_ini=nodes_Tair;
866nodes_q_ini=nodes_q;
867}
868
869#pragma GCC diagnostic pop

Archive Download this file

Revision: HEAD