GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 80.9% 233 / 0 / 288
Functions: 89.5% 17 / 0 / 19
Branches: 66.0% 103 / 0 / 156

snow.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // snow.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 11/01/21 (Build 5.2.0)
7 // Author: L. Rossman
8 //
9 // Models snow melt processes.
10 //
11 // Update History
12 // ==============
13 // Build 5.1.008:
14 // - Adjustment of snowmelt and subcatchment's net precipitation for area
15 // covered by snow was corrected.
16 // - Area covered by snow now included in calculation of rate that liquid
17 // water leaves a snowpack.
18 // Build 5.2.0:
19 // - Subcatchment snow pack area should not include LID area.
20 //-----------------------------------------------------------------------------
21 #define _CRT_SECURE_NO_DEPRECATE
22
23 #include <stdlib.h>
24 #include <string.h>
25 #include <math.h>
26 #include "headers.h"
27
28 //-----------------------------------------------------------------------------
29 // Constants
30 //-----------------------------------------------------------------------------
31 // These symbolize the keywords listed in SnowmeltWords in keywords.c
32 enum SnowKeywords {SNOW_PLOWABLE, SNOW_IMPERV, SNOW_PERV, SNOW_REMOVAL};
33
34 //-----------------------------------------------------------------------------
35 // External functions (declared in funcs.h)
36 //-----------------------------------------------------------------------------
37 // snow_createSnowpack (called from subcatch_setParams)
38 // snow_initSnowpack (called from subcatch_initState)
39 // snow_initSnowmelt (called from createObjects in project.c)
40 // snow_validateSnowmelt(called from project_validate)
41 // snow_readMeltParams (called from parseLine in input.c)
42 // snow_setMeltCoeffs (called from setTemp in climate.c)
43 // snow_plowSnow (called from runoff_execute)
44 // snow_getSnowMelt (called from subcatch_getRunoff)
45 // snow_getSnowCover (called from massbal_open)
46 // snow_getState (called from saveRunoff in hotstart.c)
47
48 //-----------------------------------------------------------------------------
49 // Local functions
50 //-----------------------------------------------------------------------------
51 static void setMeltParams(int i, int k, double x[]);
52 static double getRainmelt(double rainfall);
53 static double getArealDepletion(TSnowpack* snowpack, int i, double snowfall,
54 double tStep);
55 static double getArealSnowCover(int i, double awesi);
56 static double meltSnowpack(TSnowpack* snowpack, int i, double rmelt, double asc,
57 double snowfall, double tStep);
58 static double reduceColdContent(TSnowpack* snowpack, int i, double smelt,
59 double ccFactor);
60 static double routeSnowmelt(TSnowpack* snowpack, int i, double smelt, double asc,
61 double rainfall, double tStep);
62 static void updateColdContent(TSnowpack* snowpack, int i, double asc,
63 double snowfall, double tStep);
64
65
66 //=============================================================================
67
68 12 int snow_readMeltParams(char* tok[], int ntoks)
69 //
70 // Input: tok[] = array of string tokens
71 // ntoks = number of tokens
72 // Output: returns error code
73 // Purpose: reads snow melt parameters from a tokenized line of input data.
74 //
75 // Format of data are:
76 // Name SubArea Cmin Cmax Tbase FWF SD0 FW0 SNN0/SD100
77 // Name REMOVAL SDplow Fout Fimperv Fperv Fimelt Fsubcatch (Subcatch)
78 //
79 {
80 int i, j, k, m, n;
81 double x[7];
82
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
12 if ( ntoks < 8 ) return error_setInpError(ERR_ITEMS, "");
83
84 // --- save snow melt parameter set name if not already done so
85 12 j = project_findObject(SNOWMELT, tok[0]);
86
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
12 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
87
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 9 times.
12 if ( Snowmelt[j].ID == NULL )
88 3 Snowmelt[j].ID = project_findID(SNOWMELT, tok[0]);
89
90 // --- identify data keyword
91 12 k = findmatch(tok[1], SnowmeltWords);
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
12 if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[1]);
93
94 // --- number of parameters to read
95 12 n = 7; // 7 for subareas
96
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 9 times.
12 if ( k == SNOW_REMOVAL ) n = 6; // 6 for Removal
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
12 if ( ntoks < n + 2 ) return error_setInpError(ERR_ITEMS, "");
98
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 12 times.
96 for (i=0; i<7; i++) x[i] = 0.0;
99
100 // --- parse each parameter
101
2/2
✓ Branch 0 taken 81 times.
✓ Branch 1 taken 12 times.
93 for (i=0; i<n; i++)
102 {
103
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 81 times.
81 if ( ! getDouble(tok[i+2], &x[i]) )
104 return error_setInpError(ERR_NUMBER, tok[i+2]);
105 }
106
107 // --- parse name of subcatch receiving snow plowed from current subcatch
108
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 9 times.
12 if ( k == SNOW_REMOVAL )
109 {
110 3 x[6] = -1.0;
111
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if ( ntoks >= 9 )
112 {
113 m = project_findObject(SUBCATCH, tok[8]);
114 if ( m < 0 ) return error_setInpError(ERR_NAME, tok[8]);
115 x[6] = m;
116 }
117 }
118
119 // --- save snow melt parameters
120 12 setMeltParams(j, k, x);
121 12 return 0;
122 }
123
124 //=============================================================================
125
126 3 int snow_createSnowpack(int j, int k)
127 //
128 // Input: j = subcatchment index
129 // k = snow melt parameter set index
130 // Output: returns TRUE if successful
131 // Purpose: creates a snowpack object for a subcacthment.
132 //
133 {
134 TSnowpack* snowpack;
135 3 snowpack = (TSnowpack *) malloc(sizeof(TSnowpack));
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if ( !snowpack ) return FALSE;
137 3 Subcatch[j].snowpack = snowpack;
138 3 snowpack->snowmeltIndex = k;
139 3 return TRUE;
140 }
141
142 //=============================================================================
143
144 3 void snow_initSnowpack(int j)
145 //
146 // Input: j = subcatchment index
147 // Output: none
148 // Purpose: initializes state of a subcatchment's snow pack.
149 //
150 {
151 int i; // snow sub-area index
152 int k; // snowmelt parameter set index
153 double f; // fraction of impervious area plowable
154 3 double snowDepth = 0.0; // snow depth on entire subcatchment (ft)
155 TSnowpack* snowpack; // ptr. to snow pack object
156
157 // --- get ptr. to subcatchment's snow pack object
158 3 snowpack = Subcatch[j].snowpack;
159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if ( snowpack == NULL ) return;
160
161 // --- identify index of snow melt data set used by snow pack
162 3 k = Subcatch[j].snowpack->snowmeltIndex;
163
164 // --- find fractional area of each snow surface
165 3 f = Snowmelt[k].snn;
166 3 snowpack->fArea[SNOW_PLOWABLE] = f * Subcatch[j].fracImperv;
167 3 snowpack->fArea[SNOW_IMPERV] = (1.0 - f) * Subcatch[j].fracImperv;
168 3 snowpack->fArea[SNOW_PERV] = 1.0 - Subcatch[j].fracImperv;
169
170 // --- initialize state of snow pack on each snow surface
171
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
12 for (i=SNOW_PLOWABLE; i<=SNOW_PERV; i++)
172 {
173
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
9 if ( snowpack->fArea[i] > 0.0 )
174 {
175 6 snowpack->wsnow[i] = Snowmelt[k].wsnow[i];
176 6 snowpack->fw[i] = Snowmelt[k].fwnow[i];
177 }
178 else
179 {
180 3 snowpack->wsnow[i] = 0.0;
181 3 snowpack->fw[i] = 0.0;
182 }
183 9 snowpack->coldc[i] = 0.0;
184 9 snowpack->ati[i] = Snowmelt[k].tbase[i];
185 9 snowpack->awe[i] = 1.0;
186 9 snowDepth += snowpack->wsnow[i] * snowpack->fArea[i];
187 }
188 3 Subcatch[j].newSnowDepth = snowDepth;
189 }
190
191 //=============================================================================
192
193 3 void snow_initSnowmelt(int j)
194 //
195 // Input: j = snowmelt parameter set index
196 // Output: none
197 // Purpose: initializes values in a snow melt parameter set.
198 //
199 {
200 int i, k;
201
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
12 for (i=0; i<3; i++)
202 {
203 9 Snowmelt[j].snn = 0.0;
204 9 Snowmelt[j].si[i] = 0.0;
205 9 Snowmelt[j].dhmin[i] = 0.0;
206 9 Snowmelt[j].dhmax[i] = 0.0;
207 9 Snowmelt[j].tbase[i] = 0.0;
208 9 Snowmelt[j].fwfrac[i] = 0.0;
209 9 Snowmelt[j].wsnow[i] = 0.0;
210 9 Snowmelt[j].fwnow[i] = 0.0;
211 9 Snowmelt[j].weplow = 1.0e6;
212
2/2
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 9 times.
54 for (k=0; k<5; k++) Snowmelt[j].sfrac[k] = 0.0;
213 9 Snowmelt[j].toSubcatch = -1;
214 }
215 3 }
216
217 //=============================================================================
218
219 3 void snow_validateSnowmelt(int j)
220 //
221 // Input: j = snowmelt parameter set index
222 // Output: none
223 // Purpose: checks for valid values in a snow melt parameter set.
224 //
225 {
226 int k;
227 3 char err = FALSE;
228 3 double sum = 0.0;
229
230
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
12 for ( k = SNOW_PLOWABLE; k <= SNOW_PERV; k++ )
231 {
232 // --- check melt coeffs.
233
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
9 if ( Snowmelt[j].dhmin[k] > Snowmelt[j].dhmax[k] ) err = TRUE;
234
235 // --- check free water fraction
236
1/2
✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
9 if ( Snowmelt[j].fwfrac[k] < 0.0 ||
237
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
9 Snowmelt[j].fwfrac[k] > 1.0) err = TRUE;
238 }
239
240 // --- check fraction of imperv. area plowable
241
2/4
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
3 if ( Snowmelt[j].snn < 0.0 || Snowmelt[j].snn > 1.0 ) err = TRUE;
242
243 // --- check that removal fractions sum <= 1.0
244
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 3 times.
18 for ( k=0; k<5; k++ ) sum += Snowmelt[j].sfrac[k];
245
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if ( sum > 1.01 ) err = TRUE;
246
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if ( err ) report_writeErrorMsg(ERR_SNOWPACK_PARAMS, Snowmelt[j].ID);
247 3 }
248
249 //=============================================================================
250
251 void snow_getState(int i, int j, double x[])
252 //
253 // Input: i = subcatchment index
254 // j = snow pack sub-area index
255 // Output: updates array of snow pack state variables x
256 // Purpose: retrieves the current state of a snow pack object.
257 //
258 {
259 TSnowpack* snowpack = Subcatch[i].snowpack;
260 if ( snowpack == NULL ) return;
261 x[0] = snowpack->wsnow[j];
262 x[1] = snowpack->fw[j];
263 x[2] = snowpack->coldc[j];
264 x[3] = snowpack->ati[j];
265 x[4] = snowpack->awe[j];
266 }
267
268 //=============================================================================
269
270 void snow_setState(int i, int j, double x[])
271 //
272 // Input: i = subcatchment index
273 // j = snow pack sub-area index
274 // x = array of snow pack state variables
275 // Output: none
276 // Purpose: sets the current state of a snow pack object.
277 //
278 {
279 TSnowpack* snowpack = Subcatch[i].snowpack;
280 if ( snowpack == NULL ) return;
281 snowpack->wsnow[j] = x[0];
282 snowpack->fw[j] = x[1];
283 snowpack->coldc[j] = x[2];
284 snowpack->ati[j] = x[3];
285 snowpack->awe[j] = x[4];
286 }
287
288 //=============================================================================
289
290 12 void setMeltParams(int j, int k, double x[])
291 //
292 // Input: j = snowmelt parameter set index
293 // k = data category index
294 // x = array of snow parameter values
295 // Output: none
296 // Purpose: assigns values to parameters in a snow melt data set.
297 //
298 {
299 int i;
300
301 // --- snow pack melt parameters
302
3/4
✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 3 times.
12 if ( k >= SNOW_PLOWABLE && k <= SNOW_PERV )
303 {
304 // --- min/max melt coeffs.
305 9 Snowmelt[j].dhmin[k] = x[0] * UCF(TEMPERATURE) / UCF(RAINFALL);
306 9 Snowmelt[j].dhmax[k] = x[1] * UCF(TEMPERATURE) / UCF(RAINFALL);
307
308 // --- base melt temp (deg F)
309 9 Snowmelt[j].tbase[k] = x[2];
310
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
9 if ( UnitSystem == SI )
311 Snowmelt[j].tbase[k] = (9./5.) * Snowmelt[j].tbase[k] + 32.0;
312
313 // --- free water fractions
314 9 Snowmelt[j].fwfrac[k] = x[3];
315
316 // --- initial snow depth & free water depth
317 9 Snowmelt[j].wsnow[k] = x[4] / UCF(RAINDEPTH);
318
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 2 times.
9 x[5] = MIN(x[5], (x[3]*x[4]));
319 9 Snowmelt[j].fwnow[k] = x[5] / UCF(RAINDEPTH);
320
321 // --- fraction of impervious area that is plowable
322
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 6 times.
9 if ( k == SNOW_PLOWABLE ) Snowmelt[j].snn = x[6];
323
324 // --- min. depth for 100% areal coverage on remaining
325 // impervious area or total pervious area
326 6 else Snowmelt[j].si[k] = x[6] / UCF(RAINDEPTH);
327 }
328
329 // --- removal parameters
330
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
3 else if ( k == SNOW_REMOVAL )
331 {
332 3 Snowmelt[j].weplow = x[0] / UCF(RAINDEPTH);
333
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 3 times.
18 for (i=0; i<=4; i++) Snowmelt[j].sfrac[i] = x[i+1];
334
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 if ( x[6] >= 0.0 ) Snowmelt[j].toSubcatch = (int)(x[6] + 0.01);
335 3 else Snowmelt[j].toSubcatch = -1;
336 }
337 12 }
338
339 //=============================================================================
340
341 6 void snow_setMeltCoeffs(int j, double s)
342 //
343 // Input: j = snowmelt parameter set index
344 // s = snow season of year
345 // Output: none
346 // Purpose: sets values of snow melt coeffs. for particular time of year.
347 //
348 {
349 int k; // snow sub-area index
350
351
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
24 for (k=SNOW_PLOWABLE; k<=SNOW_PERV; k++)
352 {
353 18 Snowmelt[j].dhm[k] = 0.5 * (Snowmelt[j].dhmax[k] * (1.0 + s)
354 18 + Snowmelt[j].dhmin[k] * (1.0 - s));
355 }
356 6 }
357
358 //=============================================================================
359
360 3540 void snow_plowSnow(int j, double tStep)
361 //
362 // Input: j = subcatchment index
363 // tStep = time step (sec)
364 // Output: none
365 // Purpose: adds new snow to subcatchment and plows it between sub-areas.
366 //
367 {
368 int i; // snow sub-area index
369 int k; // snowmelt parameter set index
370 int m; // subcatchment index
371 double rainfall; // rainfall (not used)
372 double snowfall; // snowfall (ft/sec)
373 double exc; // excess snow depth (ft)
374 double f; // area ratio
375 double sfracTotal; // total fraction of snow moved
376 TSnowpack* snowpack; // ptr. to snow pack object
377
378 3540 snowpack = Subcatch[j].snowpack;
379
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3540 times.
3540 if ( !snowpack ) return;
380
381 // --- see if there's any snowfall
382 3540 gage_getPrecip(Subcatch[j].gage, &rainfall, &snowfall);
383
384 // --- add snowfall to snow pack
385
2/2
✓ Branch 0 taken 10620 times.
✓ Branch 1 taken 3540 times.
14160 for (i=SNOW_PLOWABLE; i<=SNOW_PERV; i++)
386 {
387
2/2
✓ Branch 0 taken 7080 times.
✓ Branch 1 taken 3540 times.
10620 if ( snowpack->fArea[i] > 0.0 )
388 {
389 7080 snowpack->wsnow[i] += snowfall * tStep;
390 7080 snowpack->imelt[i] = 0.0;
391 }
392 }
393
394 // --- see if there is excess snow on plowable area to remove
395
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3540 times.
3540 if ( snowpack->fArea[SNOW_PLOWABLE] > 0.0 )
396 {
397 k = snowpack->snowmeltIndex;
398 if ( snowpack->wsnow[SNOW_PLOWABLE] >= Snowmelt[k].weplow )
399 {
400 // --- excess snow to be reomoved
401 exc = snowpack->wsnow[SNOW_PLOWABLE];
402
403 // --- plow out of system
404 f = snowpack->fArea[SNOW_PLOWABLE] *
405 (Subcatch[j].area - Subcatch[j].lidArea);
406 Snow.removed += Snowmelt[k].sfrac[0] * exc * f;
407 sfracTotal = Snowmelt[k].sfrac[0];
408
409 // --- plow onto non-plowable impervious area
410 if ( snowpack->fArea[SNOW_IMPERV] > 0.0 )
411 {
412 f = snowpack->fArea[SNOW_PLOWABLE] /
413 snowpack->fArea[SNOW_IMPERV];
414 snowpack->wsnow[SNOW_IMPERV] += Snowmelt[k].sfrac[1] * exc * f;
415 sfracTotal += Snowmelt[k].sfrac[1];
416 }
417
418 // --- plow onto pervious area
419 if ( snowpack->fArea[SNOW_PERV] > 0.0 )
420 {
421 f = snowpack->fArea[SNOW_PLOWABLE] /
422 snowpack->fArea[SNOW_PERV];
423 snowpack->wsnow[SNOW_PERV] += Snowmelt[k].sfrac[2] * exc * f;
424 sfracTotal += Snowmelt[k].sfrac[2];
425 }
426
427 // --- convert to immediate melt
428 snowpack->imelt[SNOW_PLOWABLE] = Snowmelt[k].sfrac[3] * exc / tStep;
429 sfracTotal += Snowmelt[k].sfrac[3];
430
431 // --- send to another subcatchment
432 if ( Snowmelt[k].sfrac[4] > 0.0 )
433 {
434 m = Snowmelt[k].toSubcatch;
435 if ( Subcatch[m].snowpack )
436 {
437 f = Subcatch[m].snowpack->fArea[SNOW_PERV];
438 }
439 else f = 0.0;
440 if ( f > 0.0 )
441 {
442 f = snowpack->fArea[SNOW_PLOWABLE] / f;
443 Subcatch[m].snowpack->wsnow[SNOW_PERV] +=
444 Snowmelt[k].sfrac[4] * exc * f;
445 sfracTotal += Snowmelt[k].sfrac[4];
446 }
447 }
448
449 // --- reduce snow depth by amount plowed
450 sfracTotal = MIN(sfracTotal, 1.0);
451 snowpack->wsnow[SNOW_PLOWABLE] = exc * (1.0 - sfracTotal);
452 }
453 }
454 }
455
456 //=============================================================================
457
458 3540 double snow_getSnowMelt(int j, double rainfall, double snowfall, double tStep,
459 double netPrecip[])
460 //
461 // Input: j = subcatchment index
462 // rainfall = rainfall (ft/sec)
463 // snowfall = snowfall (ft/sec)
464 // tStep = time step (sec)
465 // Output: netPrecip = rainfall + snowmelt on each runoff sub-area (ft/sec),
466 // returns new snow depth over subcatchment
467 // Purpose: modifies rainfall input to subcatchment's sub-areas based on
468 // possible snow melt and updates snow depth over entire subcatchment.
469 //
470 {
471 int i; // snow sub-area index
472 double rmelt; // melt rate when rain falling (ft/sec)
473 double smelt; // snow melt from sub-area (ft/sec)
474 double asc; // frac. of sub-area snow covered
475 3540 double snowDepth = 0.0; // snow depth on entire subcatchment (ft)
476 double impervPrecip; // net precip. on imperv. area (ft/sec)
477 TSnowpack* snowpack; // ptr. to snow pack object
478
479 // --- get ptr. to subcatchment's snowpack
480 3540 snowpack = Subcatch[j].snowpack;
481
482 // --- compute snowmelt over entire subcatchment when rain falling
483 3540 rmelt = getRainmelt(rainfall);
484
485 // --- compute snow melt from each type of subarea
486
2/2
✓ Branch 0 taken 10620 times.
✓ Branch 1 taken 3540 times.
14160 for (i=SNOW_PLOWABLE; i<=SNOW_PERV; i++)
487 {
488 // --- completely melt pack if its depth is < 0.001 inch
489
2/2
✓ Branch 0 taken 6660 times.
✓ Branch 1 taken 3960 times.
10620 if ( snowpack->wsnow[i] <= 0.001 / 12.0 )
490 {
491 6660 asc = 0.0;
492 6660 smelt = 0.0;
493 6660 snowpack->imelt[i] += (snowpack->wsnow[i] + snowpack->fw[i])
494 6660 / tStep;
495 6660 snowpack->wsnow[i] = 0.0;
496 6660 snowpack->fw[i] = 0.0;
497 6660 snowpack->coldc[i] = 0.0;
498 }
499
500 // --- otherwise compute areal depletion, find snow melt
501 // and route it through pack
502 else
503 {
504 3960 asc = getArealDepletion(snowpack, i, snowfall, tStep);
505 3960 smelt = meltSnowpack(snowpack, i, rmelt, asc, snowfall, tStep);
506 3960 smelt = routeSnowmelt(snowpack, i, smelt, asc, rainfall, tStep);
507 }
508
509 // --- find net precip. over entire subcatch area
510 10620 netPrecip[i] = smelt + snowpack->imelt[i] // snow pack melt
511 10620 + rainfall*(1.0 - asc); // rainfall on non-snow area
512
513 // --- add to total snow depth on subcatchment
514 10620 snowDepth += snowpack->wsnow[i] * snowpack->fArea[i];
515 }
516
517 // --- combine netPrecip on plowable & non-plowable imperv. areas
518
1/2
✓ Branch 0 taken 3540 times.
✗ Branch 1 not taken.
3540 if ( Subcatch[j].fracImperv > 0.0 )
519 {
520 3540 impervPrecip =
521 3540 (netPrecip[SNOW_PLOWABLE] * snowpack->fArea[SNOW_PLOWABLE] +
522 3540 netPrecip[SNOW_IMPERV] * snowpack->fArea[SNOW_IMPERV]) /
523 3540 Subcatch[j].fracImperv;
524 3540 netPrecip[IMPERV0] = impervPrecip;
525 3540 netPrecip[IMPERV1] = impervPrecip;
526 }
527 3540 return snowDepth;
528 }
529
530 //=============================================================================
531
532 192 double snow_getSnowCover(int j)
533 //
534 // Input: j = subcatchment index
535 // Output: returns volume of snow cover (ft3)
536 // Purpose: computes volume of snow on a subcatchment.
537 //
538 {
539 int i;
540 192 double snowCover = 0.0; // snow cover volume (ft3)
541 TSnowpack* snowpack; // ptr. to snowpack object
542
543 192 snowpack = Subcatch[j].snowpack;
544
2/2
✓ Branch 0 taken 186 times.
✓ Branch 1 taken 6 times.
192 if ( !snowpack ) return 0.0;
545
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
24 for (i=SNOW_PLOWABLE; i<=SNOW_PERV; i++)
546 {
547 18 snowCover += (snowpack->wsnow[i] + snowpack->fw[i]) *
548 18 snowpack->fArea[i];
549 }
550 6 return snowCover * (Subcatch[j].area - Subcatch[j].lidArea);
551 }
552
553 //=============================================================================
554
555 3960 double getArealDepletion(TSnowpack* snowpack, int i, double snowfall, double tStep)
556 //
557 // Input: snowpack = ptr. to snow pack object
558 // i = snow sub-area index
559 // snowfall = snow fall rate (ft/sec)
560 // tStep = time step (sec)
561 // Output: returns fraction of sub-area with snow cover
562 // Purpose: depletes snow covered area as snow pack melts.
563 //
564 {
565 int k; // index of snow melt parameter set
566 double asc; // fraction of area with 100% cover
567 double si; // snow depth at 100% cover
568 double awesi; // depth relative to depth at 100% cover
569 double awe;
570 double sba;
571 double sbws;
572
573 // --- plowable sub-area not subject to areal depletion
574
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3960 times.
3960 if ( i == SNOW_PLOWABLE ) return 1.0;
575 3960 k = snowpack->snowmeltIndex;
576 3960 si = Snowmelt[k].si[i];
577
578 // --- no depletion if depth zero or above SI
579
3/4
✓ Branch 0 taken 3960 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 831 times.
✓ Branch 3 taken 3129 times.
3960 if ( si == 0.0 || snowpack->wsnow[i] >= si )
580 {
581 831 snowpack->awe[i] = 1.0;
582 831 return 1.0;
583 }
584
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3129 times.
3129 if ( snowpack->wsnow[i] == 0.0 )
585 {
586 snowpack->awe[i] = 1.0;
587 return 0.0;
588 }
589
590 // --- case of new snowfall
591
2/2
✓ Branch 0 taken 769 times.
✓ Branch 1 taken 2360 times.
3129 if ( snowfall > 0.0 )
592 {
593 769 awe = (snowpack->wsnow[i] - snowfall*tStep) / si;
594
1/2
✓ Branch 0 taken 769 times.
✗ Branch 1 not taken.
769 awe = MAX(awe, 0.0);
595 769 sba = getArealSnowCover(i, awe);
596 769 sbws = awe + (0.75*snowfall*tStep) / si;
597
1/2
✓ Branch 0 taken 769 times.
✗ Branch 1 not taken.
769 sbws = MIN(sbws, 1.0);
598 769 snowpack->awe[i] = awe;
599 769 snowpack->sba[i] = sba;
600 769 snowpack->sbws[i] = sbws;
601 769 return 1.0;
602 }
603
604 // --- case of no new snow
605 else
606 {
607 2360 awe = snowpack->awe[i];
608 2360 sba = snowpack->sba[i];
609 2360 sbws = snowpack->sbws[i];
610 2360 awesi = snowpack->wsnow[i] / si;
611
612 // --- relative snow depth is below start of new snow ADC
613
2/2
✓ Branch 0 taken 1916 times.
✓ Branch 1 taken 444 times.
2360 if ( awesi < snowpack->awe[i] )
614 {
615 1916 snowpack->awe[i] = 1.0;
616 1916 asc = getArealSnowCover(i, awesi);
617 }
618
619 // --- relative snow depth is above end of new snow ADC
620
2/2
✓ Branch 0 taken 432 times.
✓ Branch 1 taken 12 times.
444 else if ( awesi >= snowpack->sbws[i] )
621 {
622 432 asc = 1.0;
623 }
624
625 // --- relative snow depth is on new snow ADC
626 else
627 {
628 12 asc = sba + (1.0 - sba) / (sbws - awe) * (awesi - awe);
629 }
630 2360 return asc;
631 }
632 }
633
634 //=============================================================================
635
636 2685 double getArealSnowCover(int i, double awesi)
637 //
638 // Input: i = snow sub-area index
639 // awesi = snow depth relative to depth at 100% snow cover
640 // Output: returns fraction of sub-area with snow cover
641 // Purpose: finds x-value on areal depletion curve (ADC) for given y-value.
642 //
643 // Note: Areal depletion curves are associated with a project's Snow
644 // data structure. They plot relative snow depth (awesi)
645 // as a function of snow covererd area fraction (asc) in 10 equal
646 // awesi increments between 0 and 1.0.
647 //
648 {
649 int k; // type of ADC (impervious or pervious)
650 int m; // interval on ADC
651 double asc1, asc2; // asc values at ends of interval
652
653 // --- determine which ADC to use
654
2/2
✓ Branch 0 taken 1342 times.
✓ Branch 1 taken 1343 times.
2685 if ( i == SNOW_IMPERV ) k = 0;
655
1/2
✓ Branch 0 taken 1343 times.
✗ Branch 1 not taken.
1343 else if ( i == SNOW_PERV ) k = 1;
656 else return 1.0;
657
658 // --- locate interval on ADC that bounds awesi
659
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2685 times.
2685 if ( awesi <= 0.0 ) return 0.0;
660
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2685 times.
2685 if ( awesi >= 0.9999 ) return 1.0;
661 2685 m = (int)(awesi*10.0 + 0.00001);
662
663 // --- get asc values for either end of interval
664 2685 asc1 = Snow.adc[k][m];
665
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2685 times.
2685 if ( m >= 9) asc2 = 1.0;
666 2685 else asc2 = Snow.adc[k][m+1];
667
668 // --- return with interpolated asc value
669 2685 return asc1 + (asc2 - asc1) / 0.1 * (awesi - 0.1*(float)m);
670 }
671
672 //=============================================================================
673
674 3960 double meltSnowpack(TSnowpack* snowpack, int i, double rmelt, double asc,
675 double snowfall, double tStep)
676 //
677 // Input: snowpack = ptr. to snow pack object
678 // i = snow sub-area index
679 // rmelt = melt rate if raining (ft/sec)
680 // asc = fraction of area covered with snow
681 // snowfall = rate of snow fall (ft/sec)
682 // tStep = time step (sec)
683 // Output: returns snow melt rate (ft/sec)
684 // Purpose: computes rate of snow melt from snow sub-area.
685 //
686 {
687 int k; // snowmelt parameter set index
688 double smelt; // melt rate over sub-area (ft/sec)
689 double ccFactor; // cold content conversion factor
690
691 // --- if raining then use result found from getRainMelt
692 3960 k = snowpack->snowmeltIndex;
693
2/2
✓ Branch 0 taken 1538 times.
✓ Branch 1 taken 2422 times.
3960 if ( rmelt > 0.0 ) smelt = rmelt;
694
695 // --- else if air temp. >= base melt temp. then use degree-day eqn.
696
2/2
✓ Branch 0 taken 938 times.
✓ Branch 1 taken 1484 times.
2422 else if ( Temp.ta >= Snowmelt[k].tbase[i] )
697 {
698 938 smelt = Snowmelt[k].dhm[i] * (Temp.ta - Snowmelt[k].tbase[i]);
699 }
700
701 // --- otherwise alter cold content and return 0
702 else
703 {
704 1484 updateColdContent(snowpack, i, asc, snowfall, tStep);
705 1484 return 0.0;
706 }
707
708 // --- adjust snowmelt for area of snow cover
709 2476 smelt *= asc;
710
711 // --- reduce cold content of melting pack
712 2476 ccFactor = tStep * Snow.rnm * asc;
713 2476 smelt = reduceColdContent(snowpack, i, smelt, ccFactor);
714 2476 snowpack->ati[i] = Snowmelt[k].tbase[i];
715 2476 return smelt;
716 }
717
718 //=============================================================================
719
720 3540 double getRainmelt(double rainfall)
721 //
722 // Input: rainfall = rainfall rate (ft/sec)
723 // Output: returns snow melt rate (ft/sec)
724 // Purpose: computes rate of snow melt when rainfall occurs.
725 //
726 {
727 double uadj; // adjusted wind speed
728 double t1, t2, t3;
729 double smelt; // snow melt in in/hr
730
731 3540 rainfall = rainfall * 43200.0; // convert rain to in/hr
732
2/2
✓ Branch 0 taken 1514 times.
✓ Branch 1 taken 2026 times.
3540 if ( rainfall > 0.02 )
733 {
734 1514 uadj = 0.006 * Wind.ws;
735 1514 t1 = Temp.ta - 32.0;
736 1514 t2 = 7.5 * Temp.gamma * uadj;
737 1514 t3 = 8.5 * uadj * (Temp.ea - 0.18);
738 1514 smelt = t1 * (0.001167 + t2 + 0.007 * rainfall) + t3;
739 1514 return smelt / 43200.0;
740 }
741 2026 else return 0.0;
742 }
743
744 //=============================================================================
745
746 1484 void updateColdContent(TSnowpack* snowpack, int i, double asc, double snowfall,
747 double tStep)
748 //
749 // Input: snowpack = ptr. to snow pack object
750 // i = snow sub-area index
751 // asc = fraction of area snow covered
752 // snowfall = snow fall rate (ft/sec)
753 // tStep = time step (sec)
754 // Output: none
755 // Purpose: updates cold content of snow pack under non-melting conditions.
756 //
757 {
758 int k; // snowmelt parameter set index
759 double ati; // antecdent temperature index (deg F)
760 double cc; // snow pack cold content (ft)
761 double ccMax; // max. possible cold content (ft)
762 double tipm; // adjusted ATI weighting factor
763
764 // --- retrieve ATI & CC from snow pack object
765 1484 ati = snowpack->ati[i];
766 1484 cc = snowpack->coldc[i];
767
768 // --- if snowing, ATI = snow (air) temperature
769
2/2
✓ Branch 0 taken 220 times.
✓ Branch 1 taken 1264 times.
1484 if ( snowfall * 43200.0 > 0.02) ati = Temp.ta;
770 else
771 {
772 // convert ATI weighting factor from 6-hr to tStep time basis
773 1264 tipm = 1.0 - pow(1.0 - Snow.tipm, tStep / (6.0*3600.0));
774
775 // update ATI
776 1264 ati += tipm * (Temp.ta - ati);
777 }
778
779 // --- ATI cannot exceed snow melt base temperature
780 1484 k = snowpack->snowmeltIndex;
781
1/2
✓ Branch 0 taken 1484 times.
✗ Branch 1 not taken.
1484 ati = MIN(ati, Snowmelt[k].tbase[i]);
782
783 // --- update cold content
784 1484 cc += Snow.rnm * Snowmelt[k].dhm[i] * (ati - Temp.ta) * tStep * asc;
785
2/2
✓ Branch 0 taken 1480 times.
✓ Branch 1 taken 4 times.
1484 cc = MAX(cc, 0.0);
786
787 // --- maximum cold content based on assumed specific heat of snow
788 // of 0.007 in. water equiv. per deg. F
789 1484 ccMax = snowpack->wsnow[i] * 0.007 / 12.0 * (Snowmelt[k].tbase[i] - ati);
790
2/2
✓ Branch 0 taken 152 times.
✓ Branch 1 taken 1332 times.
1484 cc = MIN(cc, ccMax);
791
792 // --- assign updated values to snowpack
793 1484 snowpack->coldc[i] = cc;
794 1484 snowpack->ati[i] = ati;
795 1484 }
796
797 //=============================================================================
798
799 2476 double reduceColdContent(TSnowpack* snowpack, int i, double smelt, double ccFactor)
800 //
801 // Input: snowpack = ptr. to snowpack object
802 // i = snow sub-area index
803 // smelt = potential melt rate (ft/sec)
804 // ccFactor = cold content conversion factor
805 // Output: returns snow melt rate (ft/sec)
806 // Purpose: reduces cold content of snow pack adjusting melt rate accordingly.
807 //
808 {
809 double cc; // cold content of snow pack (ft)
810
811 2476 cc = snowpack->coldc[i];
812
2/2
✓ Branch 0 taken 2458 times.
✓ Branch 1 taken 18 times.
2476 if ( smelt * ccFactor > cc )
813 {
814 2458 smelt -= cc / ccFactor;
815 2458 cc = 0.0;
816 }
817 else
818 {
819 18 cc -= smelt * ccFactor;
820 18 smelt = 0.0;
821 }
822 2476 snowpack->coldc[i] = cc;
823 2476 return smelt;
824 }
825
826 //=============================================================================
827
828 3960 double routeSnowmelt(TSnowpack* snowpack, int i, double smelt, double asc,
829 double rainfall, double tStep)
830 //
831 // Input: snowpack = ptr. to snowpack object
832 // i = snow sub-area index
833 // smelt = snow melt rate (ft/sec)
834 // asc = fraction of area snow covered
835 // rainfall = rainfall rate (ft/sec)
836 // tStep = time step (sec)
837 // Output: returns rate of liquid snow melt leaving a snow pack (ft/sec)
838 // Purpose: routes snow melt through free water holding capacity of snow pack.
839 //
840 {
841 int k; // snowmelt parameter set index
842 double vmelt; // snow melt volume (ft)
843
844 // --- get volume of snowmelt over time step
845 3960 k = snowpack->snowmeltIndex;
846 3960 vmelt = smelt * tStep;
847
1/2
✓ Branch 0 taken 3960 times.
✗ Branch 1 not taken.
3960 vmelt = MIN(vmelt, snowpack->wsnow[i]);
848
849 // --- reduce snow depth by volume of snowmelt
850 3960 snowpack->wsnow[i] -= vmelt;
851
852 // --- add snowmelt volume and any rainfall on snow
853 // covered area of sub-area to snow pack's free water content
854 3960 snowpack->fw[i] += vmelt + rainfall * tStep * asc;
855
856 // --- excess free water becomes liquid melt that leaves the pack
857 3960 vmelt = snowpack->fw[i] - Snowmelt[k].fwfrac[i] * snowpack->wsnow[i];
858
2/2
✓ Branch 0 taken 2613 times.
✓ Branch 1 taken 1347 times.
3960 vmelt = MAX(vmelt, 0.0);
859
860 // --- reduce free water by liquid melt volume and return liquid melt rate
861 3960 snowpack->fw[i] -= vmelt;
862 3960 return vmelt / tStep;
863 }
864
865 //=============================================================================
866