Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/Saltation.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 contains the saltation model of Judith |
20 +------------------------------------------------------------------------------------------*/
21/********************************************************************************************/
22/*------------------------------------------------------------------------------------------+
23 | 24.08.2000: The new very complicated model. Judith says that it is my fault, if it is |
24 | wrong. I hope not ..... |
25 | Make a separate routine to model saltation only. |
26 | 26.04.2001. |
27 | Finally, on the Friday evening before the Swiss Bike Masters event, where Michael was |
28 | supposed to start, he also started to implement the last version of Judith's saltation |
29 | model. The GRID man Tuan Anh (Nguyen) had arrived and smiled. |
30 | The Gaudergrat experiment GAUDEX was in good shape and almost everything was up and |
31 | running. |
32 | 18.07.2003. |
33 +------------------------------------------------------------------------------------------*/
34/********************************************************************************************/
35
36/*------------------------------------------------------------------------------------------+
37 | Includes and Defines |
38 +------------------------------------------------------------------------------------------*/
39#include <alpine3d/snowdrift/SnowDrift.h>
40#include <alpine3d/snowdrift/checksum.h>
41#include <snowpack/libsnowpack.h>
42
43#include <math.h>
44
45using namespace std;
46
47/**
48 * @brief Calculate Saltation Fluxes for all Bottom Elements
49 * @param setbound bool
50 */
51void SnowDriftA3D::compSaltation(bool setbound)
52{
53
54unsigned int ix,iy;
55double tauS, tau_th, Ubar2;
56double weight, binding, sigma=300.;
57double dg;
58double flux_mean, cs_mean; /* What we finally want */
59
60const mio::Grid2DObject store( snowpack->getGrid(SnGrids::STORE) );
61const mio::Grid2DObject swe( snowpack->getGrid(SnGrids::SWE) );
62
63/* Calculate the Fluxes for all Bottom Elements */
64for (iy=1; iy<ny-1; iy++) {
65for (ix=1; ix<nx-1; ix++){
66if (cH.grid2D(ix,iy)==mio::IOUtils::nodata) {
67saltation(ix,iy) = 0.0; c_salt(ix,iy) = 0.0;
68continue;
69}
70//Michi: Please CHECK these two lines!
71flux_mean=0;
72cs_mean=0;
73
74if ( cH.grid2D(ix, iy) <= 0.05 || (store(ix,iy)+swe(ix,iy))<5.) {
75saltation(ix,iy) = 0.0; c_salt(ix,iy) = 0.0;
76continue;
77}
78
79/* Determine shear stress from the wind */
80/* tauS = DENSITY_AIR*DSQR(nodes[n].Km / nodes[n].lm); */ /* First try - Nice try */
81
82if (TKE) {
83/* From TKE using similarity */
84/* Attention FUDGE */
85tauS = std::min(Constants::density_air*nodes_e.grid3D(ix,iy,1)/5.,0.2);
86
87} else {
88/* From wind speed using log-profile */
89 Ubar2 = mio::Optim::pow2(nodes_u.grid3D(ix,iy,2)) + mio::Optim::pow2(nodes_v.grid3D(ix,iy,2)) + mio::Optim::pow2(nodes_w.grid3D(ix,iy,2));
90 tauS = Constants::density_air*Ubar2*mio::Optim::pow2(Saltation::karman/log((nodes_z.grid3D(ix,iy,2) - nodes_z.grid3D(ix,iy,1))/Saltation::z0_salt));
91 if (tauS == 0.) {
92printf("\n tauS zero: ix:%d, iy:%d, u:%f, v:%f"
93,ix,iy,nodes_u.grid3D(ix,iy,2),nodes_v.grid3D(ix,iy,2));
94 }
95}
96if (thresh_snow) {
97weight = 0.02*Constants::density_ice*(sp.grid2D(ix, iy) + 1.)*mio::Cst::gravity*MM_TO_M(rg.grid2D(ix, iy));
98binding = 0.0015*sigma*N3.grid2D(ix, iy)*rb.grid2D(ix, iy)*rb.grid2D(ix, iy)/rg.grid2D(ix, iy)/rg.grid2D(ix, iy);
99tau_th = std::max(tau_thresh, SnowDrift::schmidt_drift_fudge*(weight + binding));
100dg = std::min(grain_size,std::max(0.3*grain_size,2.*rg.grid2D(ix, iy)));
101} else {
102tau_th = tau_thresh; /* Pa */
103dg = grain_size;
104}/* m */
105
106/* Calculate Flux and Lower Concentration Boundary Condition*/
107if (!saltation_obj.compSaltation(tauS, tau_th, nodes_slope.grid3D(ix,iy,1)*(180./Constants::pi), dg, flux_mean, cs_mean)) {
108cout<<" Could not calculate Saltation"<<endl;
109 return;
110}
111saltation(ix,iy) = flux_mean; c_salt(ix,iy) = c_red*cs_mean;
112} /* for ix */
113}
114
115/* First set Zero Gradient Boundary Condition */
116if (setbound) {
117for (iy=0; iy<ny; iy++) {
118saltation(0,iy) = saltation(1,iy); saltation(nx-1,iy) = saltation(nx-2,iy);
119c_salt(0,iy) = c_salt(1,iy); c_salt(nx-1,iy) = c_salt(nx-2,iy);
120}
121
122//Michi: the old stuffs do not set the boundary correctly!
123for (ix=0; ix<nx; ix++) {
124saltation(ix,0) = saltation(ix,1); saltation(ix,ny-1) = saltation(ix,ny-2);
125c_salt(ix,0) = c_salt(ix,1); c_salt(ix,ny-1) = c_salt(ix,ny-2);
126}
127}
128/* End of Saltation Flux */
129// DEBUG("Checksum saltation=%lf, c_salt=%lf",checksum(saltation),checksum(c_salt));
130return;
131}
132

Archive Download this file

Revision: HEAD