Skip to content

Commit 63297df

Browse files
HishamEldardiryHishamEldardiry
authored
Add files via upload
1 parent d9cab01 commit 63297df

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

52 files changed

+29092
-0
lines changed

‎Rainfall_Runoff/SnowPackEnergyBalance.c‎

Lines changed: 329 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,329 @@
1+
/*
2+
* SUMMARY: SnowPackEnergyBalance.c - Calculate snow pack energy balance
3+
* USAGE: Part of DHSVM
4+
*
5+
* AUTHOR: Bart Nijssen
6+
* ORG: University of Washington, Department of Civil Engineering
7+
* E-MAIL: [email protected]
8+
* ORIG-DATE: 8-Oct-1996 at 09:09:29
9+
* LAST-MOD: Mon Apr 21 16:39:04 2003 by Keith Cherkauer <[email protected]>
10+
* DESCRIPTION: Calculate snow pack energy balance
11+
* DESCRIP-END.
12+
* FUNCTIONS: SnowPackEnergyBalance()
13+
* COMMENTS:
14+
*/
15+
16+
#include <stdarg.h>
17+
#include <stdio.h>
18+
#include <stdlib.h>
19+
#include <vicNl.h>
20+
21+
static char vcid[] = "$Id$";
22+
23+
/*****************************************************************************
24+
Function name: SnowPackEnergyBalance()
25+
26+
Purpose : Calculate the surface energy balance for the snow pack
27+
28+
Required :
29+
double TSurf - new estimate of effective surface temperature
30+
va_list ap - Argument list initialized by va_start(). For
31+
elements of list and order, see beginning of
32+
routine
33+
34+
Returns :
35+
double RestTerm - Rest term in the energy balance
36+
37+
Modifies :
38+
double *RefreezeEnergy - Refreeze energy (W/m2)
39+
double *vapor_flux - Mass flux of water vapor to or from the
40+
intercepted snow (m/timestep)
41+
double *blowing_flux - Mass flux of water vapor from blowing snow (m/timestep)
42+
double *surface_flux - Mass flux of water vapor from pack snow (m/timestep)
43+
44+
Comments :
45+
Reference: Bras, R. A., Hydrology, an introduction to hydrologic
46+
science, Addisson Wesley, Inc., Reading, etc., 1990.
47+
48+
Modifications:
49+
10-6-2000 modified to handle partial snow cover including the
50+
advection of sensible heat from bare patches to the edges
51+
of the remaining snow cover. KAC
52+
11-18-02 Modified to include the effects of blowing snow on the
53+
accumulation and ablation of the snowpack. LCB
54+
04-21-03 Removed constant variable declarations. Those still
55+
required by the VIC model are now defined in
56+
vicNl_def.h. KAC
57+
16-Jul-04 Renamed VaporMassFlux, BlowingMassFlux, and SurfaceMassFlux
58+
to vapor_flux, blowing_flux, and surface_flux, respectively,
59+
to denote fact that their units are m/timestep rather than
60+
kg/m2s. Created new variables VaporMassFlux, BlowingMassFlux,
61+
and SurfaceMassFlux with units of kg/m2s. The addresses of
62+
the *MassFlux variables are passed to latent_heat_from_snow()
63+
where values for the variables are computed. After these
64+
values are computed, vapor_flux, blowing_flux and surface_flux
65+
are derived from them by unit conversion. vapor_flux,
66+
blowing_flux, and surface_flux are the variables that are
67+
passed in/out of this function. TJB
68+
16-Jul-04 Changed the type of the last few variables (lag_one, iveg,
69+
etc) in the va_list to be double. For some reason, passing
70+
them as float or int caused them to become garbage. This may
71+
have to do with the fact that they followed variables of type
72+
(double *) in va_list, which may have caused memory alignment
73+
problems. TJB
74+
05-Aug-04 Removed lag_one, sigma_slope, fetch, iveg, Nveg, month, overstory,
75+
LastSnow, and SnowDepth from argument list, since they were only
76+
needed to pass to latent_heat_from_snow(), which no longer needs
77+
them. TJB
78+
28-Sep-04 Added Ra_used to store the aerodynamic resistance used in flux
79+
calculations. TJB
80+
2006-Sep-23 Replaced redundant STEFAN_B constant with STEFAN_B_B. TJB
81+
2009-Jan-16 Modified aero_resist_used and Ra_used to become arrays of
82+
two elements (surface and overstory); added
83+
options.AERO_RESIST_CANSNOW. TJB
84+
2009-Sep-19 Added Added ground flux computation consistent with 4.0.6. TJB
85+
2013-Dec-27 Moved SPATIAL_SNOW from compile-time to run-time options. TJB
86+
87+
*****************************************************************************/
88+
double SnowPackEnergyBalance(double TSurf, va_list ap)
89+
{
90+
91+
extern option_struct options;
92+
93+
const char *Routine = "SnowPackEnergyBalance";
94+
95+
/* Define Variable Argument List */
96+
97+
/* General Model Parameters */
98+
double Dt; /* Model time step (sec) */
99+
double Ra; /* Aerodynamic resistance (s/m) */
100+
double *Ra_used; /* Aerodynamic resistance (s/m) after stability correction */
101+
102+
/* Vegetation Parameters */
103+
double Displacement; /* Displacement height (m) */
104+
double Z; /* Reference height (m) */
105+
double *Z0; /* surface roughness height (m) */
106+
107+
/* Atmospheric Forcing Variables */
108+
double AirDens; /* Density of air (kg/m3) */
109+
double EactAir; /* Actual vapor pressure of air (Pa) */
110+
double LongSnowIn; /* Incoming longwave radiation (W/m2) */
111+
double Lv; /* Latent heat of vaporization (J/kg3) */
112+
double Press; /* Air pressure (Pa) */
113+
double Rain; /* Rain fall (m/timestep) */
114+
double NetShortUnder; /* Net incident shortwave radiation
115+
(W/m2) */
116+
double Vpd; /* Vapor pressure deficit (Pa) */
117+
double Wind; /* Wind speed (m/s) */
118+
119+
/* Snowpack Variables */
120+
double OldTSurf; /* Surface temperature during previous time
121+
step */
122+
double SnowCoverFract; /* Fraction of area covered by snow */
123+
double SnowDepth; /* Depth of snowpack (m) */
124+
double SnowDensity; /* Density of snowpack (kg/m^3) */
125+
double SurfaceLiquidWater; /* Liquid water in the surface layer (m) */
126+
double SweSurfaceLayer; /* Snow water equivalent in surface layer
127+
(m) */
128+
129+
/* Energy Balance Components */
130+
double Tair; /* Canopy air / Air temperature (C) */
131+
double TGrnd; /* Ground surface temperature (C) */
132+
133+
double *AdvectedEnergy; /* Energy advected by precipitation (W/m2) */
134+
double *AdvectedSensibleHeat; /* Sensible heat advected from snow-free
135+
area into snow covered area (W/m^2) */
136+
double *DeltaColdContent; /* Change in cold content of surface
137+
layer (W/m2) */
138+
double *GroundFlux; /* Ground Heat Flux (W/m2) */
139+
double *LatentHeat; /* Latent heat exchange at surface (W/m2) */
140+
double *LatentHeatSub; /* Latent heat of sublimation exchange at
141+
surface (W/m2) */
142+
double *NetLongUnder; /* Net longwave radiation at snowpack
143+
surface (W/m^2) */
144+
double *RefreezeEnergy; /* Refreeze energy (W/m2) */
145+
double *SensibleHeat; /* Sensible heat exchange at surface
146+
(W/m2) */
147+
double *vapor_flux; /* Mass flux of water vapor to or from the
148+
intercepted snow (m/timestep) */
149+
double *blowing_flux; /* Mass flux of water vapor from blowing snow. (m/timestep) */
150+
double *surface_flux; /* Mass flux of water vapor from pack snow. (m/timestep) */
151+
152+
/* Internal Routine Variables */
153+
154+
double Density; /* Density of water/ice at TMean (kg/m3) */
155+
/* double LongRadOut; */ /* long wave radiation emitted by surface
156+
(W/m2) */
157+
double NetRad; /* Net radiation exchange at surface
158+
(W/m2) */
159+
double RestTerm; /* Rest term in surface energy balance
160+
(W/m2) */
161+
double TMean; /* Average temperature for time step (C) */
162+
double Tmp;
163+
double VaporMassFlux; /* Mass flux of water vapor to or from the
164+
intercepted snow (kg/m2s) */
165+
double BlowingMassFlux; /* Mass flux of water vapor from blowing snow. (kg/m2s) */
166+
double SurfaceMassFlux; /* Mass flux of water vapor from pack snow. (kg/m2s) */
167+
168+
/* Assign the elements of the array to the appropriate variables. The list
169+
is traversed as if the elements are doubles, because:
170+
171+
In the variable-length part of variable-length argument lists, the old
172+
``default argument promotions'' apply: arguments of type double are
173+
always promoted (widened) to type double, and types char and short int
174+
are promoted to int. Therefore, it is never correct to invoke
175+
va_arg(argp, double); instead you should always use va_arg(argp,
176+
double).
177+
178+
(quoted from the comp.lang.c FAQ list)
179+
*/
180+
181+
/* General Model Parameters */
182+
Dt = (double) va_arg(ap, double);
183+
Ra = (double) va_arg(ap, double);
184+
Ra_used = (double *) va_arg(ap, double *);
185+
186+
/* Vegetation Parameters */
187+
Displacement = (double) va_arg(ap, double);
188+
Z = (double) va_arg(ap, double);
189+
Z0 = (double *) va_arg(ap, double *);
190+
191+
/* Atmospheric Forcing Variables */
192+
AirDens = (double) va_arg(ap, double);
193+
EactAir = (double) va_arg(ap, double);
194+
LongSnowIn = (double) va_arg(ap, double);
195+
Lv = (double) va_arg(ap, double);
196+
Press = (double) va_arg(ap, double);
197+
Rain = (double) va_arg(ap, double);
198+
NetShortUnder = (double) va_arg(ap, double);
199+
Vpd = (double) va_arg(ap, double);
200+
Wind = (double) va_arg(ap, double);
201+
202+
/* Snowpack Variables */
203+
OldTSurf = (double) va_arg(ap, double);
204+
SnowCoverFract = (double) va_arg(ap, double);
205+
SnowDepth = (double) va_arg(ap, double);
206+
SnowDensity = (double) va_arg(ap, double);
207+
SurfaceLiquidWater = (double) va_arg(ap, double);
208+
SweSurfaceLayer = (double) va_arg(ap, double);
209+
210+
/* Energy Balance Components */
211+
Tair = (double) va_arg(ap, double);
212+
TGrnd = (double) va_arg(ap, double);
213+
214+
AdvectedEnergy = (double *) va_arg(ap, double *);
215+
AdvectedSensibleHeat = (double *)va_arg(ap, double *);
216+
DeltaColdContent = (double *) va_arg(ap, double *);
217+
GroundFlux = (double *) va_arg(ap, double *);
218+
LatentHeat = (double *) va_arg(ap, double *);
219+
LatentHeatSub = (double *) va_arg(ap, double *);
220+
NetLongUnder = (double *) va_arg(ap, double *);
221+
RefreezeEnergy = (double *) va_arg(ap, double *);
222+
SensibleHeat = (double *) va_arg(ap, double *);
223+
vapor_flux = (double *) va_arg(ap, double *);
224+
blowing_flux = (double *) va_arg(ap, double *);
225+
surface_flux = (double *) va_arg(ap, double *);
226+
227+
/* Calculate active temp for energy balance as average of old and new */
228+
229+
TMean = TSurf;
230+
Density = RHO_W;
231+
232+
/* Correct aerodynamic conductance for stable conditions
233+
Note: If air temp >> snow temp then aero_cond -> 0 (i.e. very stable)
234+
velocity (vel_2m) is expected to be in m/sec */
235+
236+
/* Apply the stability correction to the aerodynamic resistance
237+
NOTE: In the old code 2m was passed instead of Z-Displacement. I (bart)
238+
think that it is more correct to calculate ALL fluxes at the same
239+
reference level */
240+
241+
242+
if (Wind > 0.0)
243+
Ra_used[0] = Ra / StabilityCorrection(Z, 0.f, TMean, Tair, Wind, Z0[2]);
244+
else
245+
Ra_used[0] = HUGE_RESIST;
246+
247+
/* Calculate longwave exchange and net radiation */
248+
249+
Tmp = TMean + KELVIN;
250+
(*NetLongUnder) = LongSnowIn - STEFAN_B * Tmp * Tmp * Tmp * Tmp;
251+
NetRad = NetShortUnder + (*NetLongUnder);
252+
253+
/* Calculate the sensible heat flux */
254+
255+
*SensibleHeat = AirDens * Cp * (Tair - TMean) / Ra_used[0];
256+
257+
if (options.SPATIAL_SNOW) {
258+
/* Add in Sensible heat flux turbulent exchange from surrounding
259+
snow free patches - if present */
260+
if ( SnowCoverFract > 0 ) {
261+
*(AdvectedSensibleHeat) = advected_sensible_heat(SnowCoverFract,
262+
AirDens, Tair, TGrnd,
263+
Ra_used[0]);
264+
}
265+
else (*AdvectedSensibleHeat) = 0;
266+
}
267+
else {
268+
(*AdvectedSensibleHeat) = 0;
269+
}
270+
271+
/* Convert sublimation terms from m/timestep to kg/m2s */
272+
VaporMassFlux = *vapor_flux * Density / Dt;
273+
BlowingMassFlux = *blowing_flux * Density / Dt;
274+
SurfaceMassFlux = *surface_flux * Density / Dt;
275+
276+
/* Calculate the mass flux of ice to or from the surface layer */
277+
278+
/* Calculate the saturated vapor pressure in the snow pack,
279+
(Equation 3.32, Bras 1990) */
280+
281+
latent_heat_from_snow(AirDens, EactAir, Lv, Press, Ra_used[0], TMean, Vpd,
282+
LatentHeat, LatentHeatSub, &VaporMassFlux, &BlowingMassFlux,
283+
&SurfaceMassFlux);
284+
285+
/* Convert sublimation terms from kg/m2s to m/timestep */
286+
*vapor_flux = VaporMassFlux * Dt / Density;
287+
*blowing_flux = BlowingMassFlux * Dt / Density;
288+
*surface_flux = SurfaceMassFlux * Dt / Density;
289+
290+
/* Calculate advected heat flux from rain
291+
Equation 7.3.12 from H.B.H. for rain falling on melting snowpack */
292+
293+
if ( TMean == 0 )
294+
*AdvectedEnergy = (CH_WATER * (Tair) * Rain) / (Dt);
295+
else
296+
*AdvectedEnergy = 0.;
297+
298+
/* Calculate change in cold content */
299+
*DeltaColdContent = CH_ICE * SweSurfaceLayer * (TSurf - OldTSurf) /
300+
(Dt);
301+
302+
/* Calculate Ground Heat Flux */
303+
if(SnowDepth>0.) {
304+
*GroundFlux = 2.9302e-6 * SnowDensity * SnowDensity
305+
* (TGrnd - TMean) / SnowDepth / (Dt);
306+
}
307+
else *GroundFlux=0;
308+
*DeltaColdContent -= *GroundFlux;
309+
310+
/* Calculate energy balance error at the snowpack surface */
311+
RestTerm = NetRad + *SensibleHeat + *LatentHeat + *LatentHeatSub
312+
+ *AdvectedEnergy + *GroundFlux - *DeltaColdContent
313+
+ *AdvectedSensibleHeat;
314+
315+
*RefreezeEnergy = (SurfaceLiquidWater * Lf * Density)/(Dt);
316+
317+
if (TSurf == 0.0 && RestTerm > -(*RefreezeEnergy)) {
318+
*RefreezeEnergy = -RestTerm; /* available energy input over cold content
319+
used to melt, i.e. Qrf is negative value
320+
(energy out of pack)*/
321+
RestTerm = 0.0;
322+
}
323+
else {
324+
RestTerm += *RefreezeEnergy; /* add this positive value to the pack */
325+
}
326+
327+
return RestTerm;
328+
}
329+

‎Rainfall_Runoff/SnowPackEnergyBalance.o‎

15.9 KB
Binary file not shown.

0 commit comments

Comments
 (0)