GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 94.1% 271 / 0 / 288
Functions: 100.0% 19 / 0 / 19
Branches: 72.4% 113 / 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 32 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 32 times.
32 if ( ntoks < 8 ) return error_setInpError(ERR_ITEMS, "");
83
84 // --- save snow melt parameter set name if not already done so
85 32 j = project_findObject(SNOWMELT, tok[0]);
86
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
32 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
87
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 24 times.
32 if ( Snowmelt[j].ID == NULL )
88 8 Snowmelt[j].ID = project_findID(SNOWMELT, tok[0]);
89
90 // --- identify data keyword
91 32 k = findmatch(tok[1], SnowmeltWords);
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
32 if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[1]);
93
94 // --- number of parameters to read
95 32 n = 7; // 7 for subareas
96
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 24 times.
32 if ( k == SNOW_REMOVAL ) n = 6; // 6 for Removal
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
32 if ( ntoks < n + 2 ) return error_setInpError(ERR_ITEMS, "");
98
2/2
✓ Branch 0 taken 224 times.
✓ Branch 1 taken 32 times.
256 for (i=0; i<7; i++) x[i] = 0.0;
99
100 // --- parse each parameter
101
2/2
✓ Branch 0 taken 216 times.
✓ Branch 1 taken 32 times.
248 for (i=0; i<n; i++)
102 {
103
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 216 times.
216 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 8 times.
✓ Branch 1 taken 24 times.
32 if ( k == SNOW_REMOVAL )
109 {
110 8 x[6] = -1.0;
111
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
8 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 32 setMeltParams(j, k, x);
121 32 return 0;
122 }
123
124 //=============================================================================
125
126 2304 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 2304 snowpack = (TSnowpack *) malloc(sizeof(TSnowpack));
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2304 times.
2304 if ( !snowpack ) return FALSE;
137 2304 Subcatch[j].snowpack = snowpack;
138 2304 snowpack->snowmeltIndex = k;
139 2304 return TRUE;
140 }
141
142 //=============================================================================
143
144 2304 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 2304 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 2304 snowpack = Subcatch[j].snowpack;
159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2304 times.
2304 if ( snowpack == NULL ) return;
160
161 // --- identify index of snow melt data set used by snow pack
162 2304 k = Subcatch[j].snowpack->snowmeltIndex;
163
164 // --- find fractional area of each snow surface
165 2304 f = Snowmelt[k].snn;
166 2304 snowpack->fArea[SNOW_PLOWABLE] = f * Subcatch[j].fracImperv;
167 2304 snowpack->fArea[SNOW_IMPERV] = (1.0 - f) * Subcatch[j].fracImperv;
168 2304 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 6912 times.
✓ Branch 1 taken 2304 times.
9216 for (i=SNOW_PLOWABLE; i<=SNOW_PERV; i++)
172 {
173
2/2
✓ Branch 0 taken 6911 times.
✓ Branch 1 taken 1 time.
6912 if ( snowpack->fArea[i] > 0.0 )
174 {
175 6911 snowpack->wsnow[i] = Snowmelt[k].wsnow[i];
176 6911 snowpack->fw[i] = Snowmelt[k].fwnow[i];
177 }
178 else
179 {
180 1 snowpack->wsnow[i] = 0.0;
181 1 snowpack->fw[i] = 0.0;
182 }
183 6912 snowpack->coldc[i] = 0.0;
184 6912 snowpack->ati[i] = Snowmelt[k].tbase[i];
185 6912 snowpack->awe[i] = 1.0;
186 6912 snowDepth += snowpack->wsnow[i] * snowpack->fArea[i];
187 }
188 2304 Subcatch[j].newSnowDepth = snowDepth;
189 }
190
191 //=============================================================================
192
193 8 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 24 times.
✓ Branch 1 taken 8 times.
32 for (i=0; i<3; i++)
202 {
203 24 Snowmelt[j].snn = 0.0;
204 24 Snowmelt[j].si[i] = 0.0;
205 24 Snowmelt[j].dhmin[i] = 0.0;
206 24 Snowmelt[j].dhmax[i] = 0.0;
207 24 Snowmelt[j].tbase[i] = 0.0;
208 24 Snowmelt[j].fwfrac[i] = 0.0;
209 24 Snowmelt[j].wsnow[i] = 0.0;
210 24 Snowmelt[j].fwnow[i] = 0.0;
211 24 Snowmelt[j].weplow = 1.0e6;
212
2/2
✓ Branch 0 taken 120 times.
✓ Branch 1 taken 24 times.
144 for (k=0; k<5; k++) Snowmelt[j].sfrac[k] = 0.0;
213 24 Snowmelt[j].toSubcatch = -1;
214 }
215 8 }
216
217 //=============================================================================
218
219 8 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 8 char err = FALSE;
228 8 double sum = 0.0;
229
230
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 8 times.
32 for ( k = SNOW_PLOWABLE; k <= SNOW_PERV; k++ )
231 {
232 // --- check melt coeffs.
233
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
24 if ( Snowmelt[j].dhmin[k] > Snowmelt[j].dhmax[k] ) err = TRUE;
234
235 // --- check free water fraction
236
1/2
✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
24 if ( Snowmelt[j].fwfrac[k] < 0.0 ||
237
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
24 Snowmelt[j].fwfrac[k] > 1.0) err = TRUE;
238 }
239
240 // --- check fraction of imperv. area plowable
241
2/4
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
8 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 40 times.
✓ Branch 1 taken 8 times.
48 for ( k=0; k<5; k++ ) sum += Snowmelt[j].sfrac[k];
245
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
8 if ( sum > 1.01 ) err = TRUE;
246
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
8 if ( err ) report_writeErrorMsg(ERR_SNOWPACK_PARAMS, Snowmelt[j].ID);
247 8 }
248
249 //=============================================================================
250
251 987 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 987 TSnowpack* snowpack = Subcatch[i].snowpack;
260
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 987 times.
987 if ( snowpack == NULL ) return;
261 987 x[0] = snowpack->wsnow[j];
262 987 x[1] = snowpack->fw[j];
263 987 x[2] = snowpack->coldc[j];
264 987 x[3] = snowpack->ati[j];
265 987 x[4] = snowpack->awe[j];
266 }
267
268 //=============================================================================
269
270 987 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 987 TSnowpack* snowpack = Subcatch[i].snowpack;
280
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 987 times.
987 if ( snowpack == NULL ) return;
281 987 snowpack->wsnow[j] = x[0];
282 987 snowpack->fw[j] = x[1];
283 987 snowpack->coldc[j] = x[2];
284 987 snowpack->ati[j] = x[3];
285 987 snowpack->awe[j] = x[4];
286 }
287
288 //=============================================================================
289
290 32 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 32 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 8 times.
32 if ( k >= SNOW_PLOWABLE && k <= SNOW_PERV )
303 {
304 // --- min/max melt coeffs.
305 24 Snowmelt[j].dhmin[k] = x[0] * UCF(TEMPERATURE) / UCF(RAINFALL);
306 24 Snowmelt[j].dhmax[k] = x[1] * UCF(TEMPERATURE) / UCF(RAINFALL);
307
308 // --- base melt temp (deg F)
309 24 Snowmelt[j].tbase[k] = x[2];
310
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
24 if ( UnitSystem == SI )
311 Snowmelt[j].tbase[k] = (9./5.) * Snowmelt[j].tbase[k] + 32.0;
312
313 // --- free water fractions
314 24 Snowmelt[j].fwfrac[k] = x[3];
315
316 // --- initial snow depth & free water depth
317 24 Snowmelt[j].wsnow[k] = x[4] / UCF(RAINDEPTH);
318
1/2
✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
24 x[5] = MIN(x[5], (x[3]*x[4]));
319 24 Snowmelt[j].fwnow[k] = x[5] / UCF(RAINDEPTH);
320
321 // --- fraction of impervious area that is plowable
322
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 16 times.
24 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 16 else Snowmelt[j].si[k] = x[6] / UCF(RAINDEPTH);
327 }
328
329 // --- removal parameters
330
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
8 else if ( k == SNOW_REMOVAL )
331 {
332 8 Snowmelt[j].weplow = x[0] / UCF(RAINDEPTH);
333
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 8 times.
48 for (i=0; i<=4; i++) Snowmelt[j].sfrac[i] = x[i+1];
334
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
8 if ( x[6] >= 0.0 ) Snowmelt[j].toSubcatch = (int)(x[6] + 0.01);
335 8 else Snowmelt[j].toSubcatch = -1;
336 }
337 32 }
338
339 //=============================================================================
340
341 40 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 120 times.
✓ Branch 1 taken 40 times.
160 for (k=SNOW_PLOWABLE; k<=SNOW_PERV; k++)
352 {
353 120 Snowmelt[j].dhm[k] = 0.5 * (Snowmelt[j].dhmax[k] * (1.0 + s)
354 120 + Snowmelt[j].dhmin[k] * (1.0 - s));
355 }
356 40 }
357
358 //=============================================================================
359
360 3115063 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 3115063 snowpack = Subcatch[j].snowpack;
379
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3115063 times.
3115063 if ( !snowpack ) return;
380
381 // --- see if there's any snowfall
382 3115063 gage_getPrecip(Subcatch[j].gage, &rainfall, &snowfall);
383
384 // --- add snowfall to snow pack
385
2/2
✓ Branch 0 taken 9345189 times.
✓ Branch 1 taken 3115063 times.
12460252 for (i=SNOW_PLOWABLE; i<=SNOW_PERV; i++)
386 {
387
2/2
✓ Branch 0 taken 9344769 times.
✓ Branch 1 taken 420 times.
9345189 if ( snowpack->fArea[i] > 0.0 )
388 {
389 9344769 snowpack->wsnow[i] += snowfall * tStep;
390 9344769 snowpack->imelt[i] = 0.0;
391 }
392 }
393
394 // --- see if there is excess snow on plowable area to remove
395
2/2
✓ Branch 0 taken 3114643 times.
✓ Branch 1 taken 420 times.
3115063 if ( snowpack->fArea[SNOW_PLOWABLE] > 0.0 )
396 {
397 3114643 k = snowpack->snowmeltIndex;
398
1/2
✓ Branch 0 taken 3114643 times.
✗ Branch 1 not taken.
3114643 if ( snowpack->wsnow[SNOW_PLOWABLE] >= Snowmelt[k].weplow )
399 {
400 // --- excess snow to be reomoved
401 3114643 exc = snowpack->wsnow[SNOW_PLOWABLE];
402
403 // --- plow out of system
404 3114643 f = snowpack->fArea[SNOW_PLOWABLE] *
405 3114643 (Subcatch[j].area - Subcatch[j].lidArea);
406 3114643 Snow.removed += Snowmelt[k].sfrac[0] * exc * f;
407 3114643 sfracTotal = Snowmelt[k].sfrac[0];
408
409 // --- plow onto non-plowable impervious area
410
1/2
✓ Branch 0 taken 3114643 times.
✗ Branch 1 not taken.
3114643 if ( snowpack->fArea[SNOW_IMPERV] > 0.0 )
411 {
412 3114643 f = snowpack->fArea[SNOW_PLOWABLE] /
413 3114643 snowpack->fArea[SNOW_IMPERV];
414 3114643 snowpack->wsnow[SNOW_IMPERV] += Snowmelt[k].sfrac[1] * exc * f;
415 3114643 sfracTotal += Snowmelt[k].sfrac[1];
416 }
417
418 // --- plow onto pervious area
419
1/2
✓ Branch 0 taken 3114643 times.
✗ Branch 1 not taken.
3114643 if ( snowpack->fArea[SNOW_PERV] > 0.0 )
420 {
421 3114643 f = snowpack->fArea[SNOW_PLOWABLE] /
422 3114643 snowpack->fArea[SNOW_PERV];
423 3114643 snowpack->wsnow[SNOW_PERV] += Snowmelt[k].sfrac[2] * exc * f;
424 3114643 sfracTotal += Snowmelt[k].sfrac[2];
425 }
426
427 // --- convert to immediate melt
428 3114643 snowpack->imelt[SNOW_PLOWABLE] = Snowmelt[k].sfrac[3] * exc / tStep;
429 3114643 sfracTotal += Snowmelt[k].sfrac[3];
430
431 // --- send to another subcatchment
432
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3114643 times.
3114643 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
1/2
✓ Branch 0 taken 3114643 times.
✗ Branch 1 not taken.
3114643 sfracTotal = MIN(sfracTotal, 1.0);
451 3114643 snowpack->wsnow[SNOW_PLOWABLE] = exc * (1.0 - sfracTotal);
452 }
453 }
454 }
455
456 //=============================================================================
457
458 3115063 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 3115063 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 3115063 snowpack = Subcatch[j].snowpack;
481
482 // --- compute snowmelt over entire subcatchment when rain falling
483 3115063 rmelt = getRainmelt(rainfall);
484
485 // --- compute snow melt from each type of subarea
486
2/2
✓ Branch 0 taken 9345189 times.
✓ Branch 1 taken 3115063 times.
12460252 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 8212769 times.
✓ Branch 1 taken 1132420 times.
9345189 if ( snowpack->wsnow[i] <= 0.001 / 12.0 )
490 {
491 8212769 asc = 0.0;
492 8212769 smelt = 0.0;
493 8212769 snowpack->imelt[i] += (snowpack->wsnow[i] + snowpack->fw[i])
494 8212769 / tStep;
495 8212769 snowpack->wsnow[i] = 0.0;
496 8212769 snowpack->fw[i] = 0.0;
497 8212769 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 1132420 asc = getArealDepletion(snowpack, i, snowfall, tStep);
505 1132420 smelt = meltSnowpack(snowpack, i, rmelt, asc, snowfall, tStep);
506 1132420 smelt = routeSnowmelt(snowpack, i, smelt, asc, rainfall, tStep);
507 }
508
509 // --- find net precip. over entire subcatch area
510 9345189 netPrecip[i] = smelt + snowpack->imelt[i] // snow pack melt
511 9345189 + rainfall*(1.0 - asc); // rainfall on non-snow area
512
513 // --- add to total snow depth on subcatchment
514 9345189 snowDepth += snowpack->wsnow[i] * snowpack->fArea[i];
515 }
516
517 // --- combine netPrecip on plowable & non-plowable imperv. areas
518
1/2
✓ Branch 0 taken 3115063 times.
✗ Branch 1 not taken.
3115063 if ( Subcatch[j].fracImperv > 0.0 )
519 {
520 3115063 impervPrecip =
521 3115063 (netPrecip[SNOW_PLOWABLE] * snowpack->fArea[SNOW_PLOWABLE] +
522 3115063 netPrecip[SNOW_IMPERV] * snowpack->fArea[SNOW_IMPERV]) /
523 3115063 Subcatch[j].fracImperv;
524 3115063 netPrecip[IMPERV0] = impervPrecip;
525 3115063 netPrecip[IMPERV1] = impervPrecip;
526 }
527 3115063 return snowDepth;
528 }
529
530 //=============================================================================
531
532 4786 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 4786 double snowCover = 0.0; // snow cover volume (ft3)
541 TSnowpack* snowpack; // ptr. to snowpack object
542
543 4786 snowpack = Subcatch[j].snowpack;
544
2/2
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 4608 times.
4786 if ( !snowpack ) return 0.0;
545
2/2
✓ Branch 0 taken 13824 times.
✓ Branch 1 taken 4608 times.
18432 for (i=SNOW_PLOWABLE; i<=SNOW_PERV; i++)
546 {
547 13824 snowCover += (snowpack->wsnow[i] + snowpack->fw[i]) *
548 13824 snowpack->fArea[i];
549 }
550 4608 return snowCover * (Subcatch[j].area - Subcatch[j].lidArea);
551 }
552
553 //=============================================================================
554
555 1132420 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 1132420 times.
1132420 if ( i == SNOW_PLOWABLE ) return 1.0;
575 1132420 k = snowpack->snowmeltIndex;
576 1132420 si = Snowmelt[k].si[i];
577
578 // --- no depletion if depth zero or above SI
579
3/4
✓ Branch 0 taken 1132420 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5404 times.
✓ Branch 3 taken 1127016 times.
1132420 if ( si == 0.0 || snowpack->wsnow[i] >= si )
580 {
581 5404 snowpack->awe[i] = 1.0;
582 5404 return 1.0;
583 }
584
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1127016 times.
1127016 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 143716 times.
✓ Branch 1 taken 983300 times.
1127016 if ( snowfall > 0.0 )
592 {
593 143716 awe = (snowpack->wsnow[i] - snowfall*tStep) / si;
594
1/2
✓ Branch 0 taken 143716 times.
✗ Branch 1 not taken.
143716 awe = MAX(awe, 0.0);
595 143716 sba = getArealSnowCover(i, awe);
596 143716 sbws = awe + (0.75*snowfall*tStep) / si;
597
1/2
✓ Branch 0 taken 143716 times.
✗ Branch 1 not taken.
143716 sbws = MIN(sbws, 1.0);
598 143716 snowpack->awe[i] = awe;
599 143716 snowpack->sba[i] = sba;
600 143716 snowpack->sbws[i] = sbws;
601 143716 return 1.0;
602 }
603
604 // --- case of no new snow
605 else
606 {
607 983300 awe = snowpack->awe[i];
608 983300 sba = snowpack->sba[i];
609 983300 sbws = snowpack->sbws[i];
610 983300 awesi = snowpack->wsnow[i] / si;
611
612 // --- relative snow depth is below start of new snow ADC
613
2/2
✓ Branch 0 taken 167940 times.
✓ Branch 1 taken 815360 times.
983300 if ( awesi < snowpack->awe[i] )
614 {
615 167940 snowpack->awe[i] = 1.0;
616 167940 asc = getArealSnowCover(i, awesi);
617 }
618
619 // --- relative snow depth is above end of new snow ADC
620
2/2
✓ Branch 0 taken 781494 times.
✓ Branch 1 taken 33866 times.
815360 else if ( awesi >= snowpack->sbws[i] )
621 {
622 781494 asc = 1.0;
623 }
624
625 // --- relative snow depth is on new snow ADC
626 else
627 {
628 33866 asc = sba + (1.0 - sba) / (sbws - awe) * (awesi - awe);
629 }
630 983300 return asc;
631 }
632 }
633
634 //=============================================================================
635
636 311656 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 136869 times.
✓ Branch 1 taken 174787 times.
311656 if ( i == SNOW_IMPERV ) k = 0;
655
1/2
✓ Branch 0 taken 174787 times.
✗ Branch 1 not taken.
174787 else if ( i == SNOW_PERV ) k = 1;
656 else return 1.0;
657
658 // --- locate interval on ADC that bounds awesi
659
2/2
✓ Branch 0 taken 3948 times.
✓ Branch 1 taken 307708 times.
311656 if ( awesi <= 0.0 ) return 0.0;
660
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 307708 times.
307708 if ( awesi >= 0.9999 ) return 1.0;
661 307708 m = (int)(awesi*10.0 + 0.00001);
662
663 // --- get asc values for either end of interval
664 307708 asc1 = Snow.adc[k][m];
665
2/2
✓ Branch 0 taken 470 times.
✓ Branch 1 taken 307238 times.
307708 if ( m >= 9) asc2 = 1.0;
666 307238 else asc2 = Snow.adc[k][m+1];
667
668 // --- return with interpolated asc value
669 307708 return asc1 + (asc2 - asc1) / 0.1 * (awesi - 0.1*(float)m);
670 }
671
672 //=============================================================================
673
674 1132420 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 1132420 k = snowpack->snowmeltIndex;
693
2/2
✓ Branch 0 taken 4032 times.
✓ Branch 1 taken 1128388 times.
1132420 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 248809 times.
✓ Branch 1 taken 879579 times.
1128388 else if ( Temp.ta >= Snowmelt[k].tbase[i] )
697 {
698 248809 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 879579 updateColdContent(snowpack, i, asc, snowfall, tStep);
705 879579 return 0.0;
706 }
707
708 // --- adjust snowmelt for area of snow cover
709 252841 smelt *= asc;
710
711 // --- reduce cold content of melting pack
712 252841 ccFactor = tStep * Snow.rnm * asc;
713 252841 smelt = reduceColdContent(snowpack, i, smelt, ccFactor);
714 252841 snowpack->ati[i] = Snowmelt[k].tbase[i];
715 252841 return smelt;
716 }
717
718 //=============================================================================
719
720 3115063 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 3115063 rainfall = rainfall * 43200.0; // convert rain to in/hr
732
2/2
✓ Branch 0 taken 176697 times.
✓ Branch 1 taken 2938366 times.
3115063 if ( rainfall > 0.02 )
733 {
734 176697 uadj = 0.006 * Wind.ws;
735 176697 t1 = Temp.ta - 32.0;
736 176697 t2 = 7.5 * Temp.gamma * uadj;
737 176697 t3 = 8.5 * uadj * (Temp.ea - 0.18);
738 176697 smelt = t1 * (0.001167 + t2 + 0.007 * rainfall) + t3;
739 176697 return smelt / 43200.0;
740 }
741 2938366 else return 0.0;
742 }
743
744 //=============================================================================
745
746 879579 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 879579 ati = snowpack->ati[i];
766 879579 cc = snowpack->coldc[i];
767
768 // --- if snowing, ATI = snow (air) temperature
769
2/2
✓ Branch 0 taken 53101 times.
✓ Branch 1 taken 826478 times.
879579 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 826478 tipm = 1.0 - pow(1.0 - Snow.tipm, tStep / (6.0*3600.0));
774
775 // update ATI
776 826478 ati += tipm * (Temp.ta - ati);
777 }
778
779 // --- ATI cannot exceed snow melt base temperature
780 879579 k = snowpack->snowmeltIndex;
781
1/2
✓ Branch 0 taken 879579 times.
✗ Branch 1 not taken.
879579 ati = MIN(ati, Snowmelt[k].tbase[i]);
782
783 // --- update cold content
784 879579 cc += Snow.rnm * Snowmelt[k].dhm[i] * (ati - Temp.ta) * tStep * asc;
785
2/2
✓ Branch 0 taken 695483 times.
✓ Branch 1 taken 184096 times.
879579 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 879579 ccMax = snowpack->wsnow[i] * 0.007 / 12.0 * (Snowmelt[k].tbase[i] - ati);
790
2/2
✓ Branch 0 taken 333203 times.
✓ Branch 1 taken 546376 times.
879579 cc = MIN(cc, ccMax);
791
792 // --- assign updated values to snowpack
793 879579 snowpack->coldc[i] = cc;
794 879579 snowpack->ati[i] = ati;
795 879579 }
796
797 //=============================================================================
798
799 252841 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 252841 cc = snowpack->coldc[i];
812
2/2
✓ Branch 0 taken 252797 times.
✓ Branch 1 taken 44 times.
252841 if ( smelt * ccFactor > cc )
813 {
814 252797 smelt -= cc / ccFactor;
815 252797 cc = 0.0;
816 }
817 else
818 {
819 44 cc -= smelt * ccFactor;
820 44 smelt = 0.0;
821 }
822 252841 snowpack->coldc[i] = cc;
823 252841 return smelt;
824 }
825
826 //=============================================================================
827
828 1132420 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 1132420 k = snowpack->snowmeltIndex;
846 1132420 vmelt = smelt * tStep;
847
2/2
✓ Branch 0 taken 1128417 times.
✓ Branch 1 taken 4003 times.
1132420 vmelt = MIN(vmelt, snowpack->wsnow[i]);
848
849 // --- reduce snow depth by volume of snowmelt
850 1132420 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 1132420 snowpack->fw[i] += vmelt + rainfall * tStep * asc;
855
856 // --- excess free water becomes liquid melt that leaves the pack
857 1132420 vmelt = snowpack->fw[i] - Snowmelt[k].fwfrac[i] * snowpack->wsnow[i];
858
2/2
✓ Branch 0 taken 236595 times.
✓ Branch 1 taken 895825 times.
1132420 vmelt = MAX(vmelt, 0.0);
859
860 // --- reduce free water by liquid melt volume and return liquid melt rate
861 1132420 snowpack->fw[i] -= vmelt;
862 1132420 return vmelt / tStep;
863 }
864
865 //=============================================================================
866