GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 89.7% 113 / 0 / 126
Functions: 100.0% 8 / 0 / 8
Branches: 76.2% 64 / 0 / 84

surfqual.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // surfqual.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 07/14/23 (Build 5.2.4)
7 // Author: L. Rossman
8 //
9 // Subcatchment water quality functions.
10 //
11 // Update History
12 // ==============
13 // Build 5.1.008:
14 // - Pollutant surface buildup and washoff functions were moved here from
15 // subcatch.c.
16 // - Support for separate accounting of LID drain flows included.
17 // Build 5.1.014:
18 // - Fixed bug in computing effective BMP removal by LIDs.
19 // Build 5.2.4:
20 // - Set low runoff flow concentrations to zero before computing runoff
21 // mass loads rather than after so that they match wet weather mass
22 // inflows reported for conveyance system nodes.
23 //-----------------------------------------------------------------------------
24 #define _CRT_SECURE_NO_DEPRECATE
25
26 #include <math.h>
27 #include <string.h>
28 #include "headers.h"
29 #include "lid.h"
30
31 //-----------------------------------------------------------------------------
32 // Imported variables
33 //-----------------------------------------------------------------------------
34 // Declared in RUNOFF.C
35 extern double* OutflowLoad; // exported pollutant mass load
36
37 // Volumes (ft3) for a subcatchment over a time step declared in SUBCATCH.C
38 extern double Vinfil; // non-LID infiltration
39 extern double Vinflow; // non-LID precip + snowmelt + runon + ponded water
40 extern double Voutflow; // non-LID runoff to subcatchment's outlet
41 extern double VlidIn; // inflow to LID units
42 extern double VlidInfil; // infiltration from LID units
43 extern double VlidOut; // surface outflow from LID units
44 extern double VlidDrain; // drain outflow from LID units
45 extern double VlidReturn; // LID outflow returned to pervious area
46
47 //-----------------------------------------------------------------------------
48 // External functions (declared in funcs.h)
49 //-----------------------------------------------------------------------------
50 // surfqual_initState (called from subcatch_initState)
51 // surfqual_getWashoff (called from runoff_execute)
52 // surfqual_getBuildup (called from runoff_execute)
53 // surfqual_sweepBuildup (called from runoff_execute)
54 // surfqual_getWtdWashoff (called from addWetWeatherInflows in routing.c)
55
56 //-----------------------------------------------------------------------------
57 // Function declarations
58 //-----------------------------------------------------------------------------
59 static void findWashoffLoads(int j, double runoff);
60 static void findPondedLoads(int j, double tStep);
61 static void findLidLoads(int j, double tStep);
62
63 //=============================================================================
64
65 2393 void surfqual_initState(int j)
66 //
67 // Input: j = subcatchment index
68 // Output: none
69 // Purpose: initializes pollutant buildup, ponded mass, and washoff.
70 //
71 {
72 int p;
73
74 // --- initialize washoff quality
75
2/2
✓ Branch 0 taken 9294 times.
✓ Branch 1 taken 2393 times.
11687 for (p = 0; p < Nobjects[POLLUT]; p++)
76 {
77 9294 Subcatch[j].oldQual[p] = 0.0;
78 9294 Subcatch[j].newQual[p] = 0.0;
79 9294 Subcatch[j].pondedQual[p] = 0.0;
80 }
81
82 // --- initialize pollutant buildup
83 2393 landuse_getInitBuildup(Subcatch[j].landFactor, Subcatch[j].initBuildup,
84 2393 Subcatch[j].area, Subcatch[j].curbLength);
85 2393 }
86
87 //=============================================================================
88
89 2728851 void surfqual_getBuildup(int j, double tStep)
90 //
91 // Input: j = subcatchment index
92 // tStep = time step (sec)
93 // Output: none
94 // Purpose: adds to pollutant buildup on subcatchment surface.
95 //
96 {
97 int i; // land use index
98 int p; // pollutant index
99 double f; // land use fraction
100 double area; // land use area (acres or hectares)
101 double curb; // land use curb length (user units)
102 double oldBuildup; // buildup at start of time step
103 double newBuildup; // buildup at end of time step
104
105 // --- consider each landuse
106
2/2
✓ Branch 0 taken 130830 times.
✓ Branch 1 taken 2728851 times.
2859681 for (i = 0; i < Nobjects[LANDUSE]; i++)
107 {
108 // --- skip landuse if not in subcatch
109 130830 f = Subcatch[j].landFactor[i].fraction;
110
2/2
✓ Branch 0 taken 49474 times.
✓ Branch 1 taken 81356 times.
130830 if ( f == 0.0 ) continue;
111
112 // --- get land area (in acres or hectares) & curb length
113 81356 area = f * Subcatch[j].area * UCF(LANDAREA);
114 81356 curb = f * Subcatch[j].curbLength;
115
116 // --- examine each pollutant
117
2/2
✓ Branch 0 taken 161883 times.
✓ Branch 1 taken 81356 times.
243239 for (p = 0; p < Nobjects[POLLUT]; p++)
118 {
119 // --- see if snow-only buildup is in effect
120
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 161883 times.
161883 if (Pollut[p].snowOnly
121 && Subcatch[j].newSnowDepth < 0.001/12.0) continue;
122
123 // --- use land use's buildup function to update buildup amount
124 161883 oldBuildup = Subcatch[j].landFactor[i].buildup[p];
125 161883 newBuildup = landuse_getBuildup(i, p, area, curb, oldBuildup,
126 tStep);
127
1/2
✓ Branch 0 taken 161883 times.
✗ Branch 1 not taken.
161883 newBuildup = MAX(newBuildup, oldBuildup);
128 161883 Subcatch[j].landFactor[i].buildup[p] = newBuildup;
129 161883 massbal_updateLoadingTotals(BUILDUP_LOAD, p,
130 (newBuildup - oldBuildup));
131 }
132 }
133 2728851 }
134
135 //=============================================================================
136
137 2848978 void surfqual_sweepBuildup(int j, DateTime aDate)
138 //
139 // Input: j = subcatchment index
140 // aDate = current date/time
141 // Output: none
142 // Purpose: reduces pollutant buildup over a subcatchment if sweeping occurs.
143 //
144 {
145 int i; // land use index
146 int p; // pollutant index
147 double oldBuildup; // buildup before sweeping (lbs or kg)
148 double newBuildup; // buildup after sweeping (lbs or kg)
149
150 // --- no sweeping if there is snow on plowable impervious area
151
2/2
✓ Branch 0 taken 2771497 times.
✓ Branch 1 taken 77481 times.
2848978 if ( Subcatch[j].snowpack != NULL &&
152
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2771497 times.
2771497 Subcatch[j].snowpack->wsnow[IMPERV0] > MIN_TOTAL_DEPTH ) return;
153
154 // --- consider each land use
155
2/2
✓ Branch 0 taken 133653 times.
✓ Branch 1 taken 2848978 times.
2982631 for (i = 0; i < Nobjects[LANDUSE]; i++)
156 {
157 // --- skip land use if not in subcatchment
158
2/2
✓ Branch 0 taken 50526 times.
✓ Branch 1 taken 83127 times.
133653 if ( Subcatch[j].landFactor[i].fraction == 0.0 ) continue;
159
160 // --- see if land use is subject to sweeping
161
1/2
✓ Branch 0 taken 83127 times.
✗ Branch 1 not taken.
83127 if ( Landuse[i].sweepInterval == 0.0 ) continue;
162
163 // --- see if sweep interval has been reached
164 if ( aDate - Subcatch[j].landFactor[i].lastSwept >=
165 Landuse[i].sweepInterval )
166 {
167 // --- update time when last swept
168 Subcatch[j].landFactor[i].lastSwept = aDate;
169
170 // --- examine each pollutant
171 for (p = 0; p < Nobjects[POLLUT]; p++)
172 {
173 // --- reduce buildup by the fraction available
174 // times the sweeping effic.
175 oldBuildup = Subcatch[j].landFactor[i].buildup[p];
176 newBuildup = oldBuildup * (1.0 - Landuse[i].sweepRemoval *
177 Landuse[i].washoffFunc[p].sweepEffic);
178 newBuildup = MIN(oldBuildup, newBuildup);
179 newBuildup = MAX(0.0, newBuildup);
180 Subcatch[j].landFactor[i].buildup[p] = newBuildup;
181
182 // --- update mass balance totals
183 massbal_updateLoadingTotals(SWEEPING_LOAD, p,
184 oldBuildup - newBuildup);
185 }
186 }
187 }
188 }
189
190 //=============================================================================
191
192 3206654 void surfqual_getWashoff(int j, double runoff, double tStep)
193 //
194 // Input: j = subcatchment index
195 // runoff = total subcatchment runoff before internal re-routing or
196 // LID controls (ft/sec)
197 // tStep = time step (sec)
198 // Output: none
199 // Purpose: computes new runoff quality for a subcatchment.
200 //
201 // Considers three pollutant generating streams that are combined together:
202 // 1. washoff of pollutant buildup as described by the project's land use
203 // washoff functions,
204 // 2. complete mix mass balance of pollutants in surface ponding on
205 // non-LID area due to runon, wet deposition, infiltration, & evaporation,
206 // 3. wet deposition and runon over LID areas.
207 //
208 {
209 int p; // pollutant index
210 int hasOutflow; // TRUE if subcatchment has outflow
211 double cOut; // final washoff concentration (mass/ft3)
212 double massLoad; // pollut. mass load (mass)
213 double vLidRain; // rainfall volume on LID area (ft3)
214 double vLidRunon; // external runon volume to LID area (ft3)
215 double vSurfOut; // surface runoff volume leaving subcatchment (ft3)
216 double vOut1; // runoff volume prior to LID treatment (ft3)
217 double vOut2; // runoff volume after LID treatment (ft3)
218 double area; // subcatchment area (ft2)
219
220 // --- return if there is no area or no pollutants
221 3206654 area = Subcatch[j].area;
222
3/4
✓ Branch 0 taken 3193846 times.
✓ Branch 1 taken 12808 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3193846 times.
3206654 if ( Nobjects[POLLUT] == 0 || area == 0.0 ) return;
223
224 // --- find contributions from washoff, runon and wet precip. to OutflowLoad
225
2/2
✓ Branch 0 taken 12616450 times.
✓ Branch 1 taken 3193846 times.
15810296 for (p = 0; p < Nobjects[POLLUT]; p++) OutflowLoad[p] = 0.0;
226 3193846 findWashoffLoads(j, runoff);
227 3193846 findPondedLoads(j, tStep);
228 3193846 findLidLoads(j, tStep);
229
230 // --- contribution from direct rainfall on LID areas
231 3193846 vLidRain = Subcatch[j].rainfall * Subcatch[j].lidArea * tStep;
232
233 // --- contribution from upstream runon onto LID areas
234 // (only if LIDs occupy full subcatchment)
235 3193846 vLidRunon = 0.0;
236
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3193846 times.
3193846 if ( area == Subcatch[j].lidArea )
237 {
238 vLidRunon = Subcatch[j].runon * area * tStep;
239 }
240
241 // --- runoff volume before LID treatment (ft3)
242 // (Voutflow, computed in subcatch_getRunoff, is subcatchment
243 // runoff volume before LID treatment)
244 3193846 vOut1 = Voutflow + vLidRain + vLidRunon;
245
246 // --- surface runoff + LID drain flow volume leaving the subcatchment
247 // (Subcatch.newRunoff, computed in subcatch_getRunoff, includes
248 // any surface runoff reduction from LID treatment)
249 3193846 vSurfOut = Subcatch[j].newRunoff * tStep;
250 3193846 vOut2 = vSurfOut + VlidDrain;
251
252 // --- determine if subcatchment outflow is below a small cutoff
253 3193846 hasOutflow = (vOut2 > MIN_RUNOFF * area * tStep);
254
255 // --- for each pollutant
256
2/2
✓ Branch 0 taken 12616450 times.
✓ Branch 1 taken 3193846 times.
15810296 for (p = 0; p < Nobjects[POLLUT]; p++)
257 {
258 // --- convert washoff load to a concentration
259 12616450 cOut = 0.0;
260
4/4
✓ Branch 0 taken 3484603 times.
✓ Branch 1 taken 9131847 times.
✓ Branch 2 taken 1353432 times.
✓ Branch 3 taken 2131171 times.
12616450 if ( vOut1 > 0.0 && hasOutflow ) cOut = OutflowLoad[p] / vOut1;
261
262 // --- assign any difference between pre- and post-LID
263 // subcatchment outflow loads to BMP removal
264
2/2
✓ Branch 0 taken 576 times.
✓ Branch 1 taken 12615874 times.
12616450 if ( Subcatch[j].lidArea > 0.0 )
265 {
266 576 massLoad = cOut * (vOut1 - vOut2) * Pollut[p].mcf;
267
2/2
✓ Branch 0 taken 97 times.
✓ Branch 1 taken 479 times.
576 if (massLoad > 0.0)
268 97 massbal_updateLoadingTotals(BMP_REMOVAL_LOAD, p, massLoad);
269 }
270
271 // --- update subcatchment's cumulative runoff load in lbs (or kg)
272 12616450 massLoad = cOut * vOut2 * Pollut[p].mcf;
273 12616450 Subcatch[j].totalLoad[p] += massLoad;
274
275 // --- update mass balance for surface runoff load routed to a
276 // conveyance system node
277 // (loads from LID drains are accounted for below since they
278 // can go to different outlets than parent subcatchment)
279
3/4
✓ Branch 0 taken 416548 times.
✓ Branch 1 taken 12199902 times.
✓ Branch 2 taken 416548 times.
✗ Branch 3 not taken.
12616450 if ( (Subcatch[j].outNode >= 0 || Subcatch[j].outSubcatch == j) )
280 {
281 12616450 massLoad = cOut * vSurfOut * Pollut[p].mcf;
282 12616450 massbal_updateLoadingTotals(RUNOFF_LOAD, p, massLoad);
283 }
284
285 // --- save new washoff concentration
286 12616450 Subcatch[j].newQual[p] = cOut / LperFT3;
287 }
288
289 // --- add contribution of LID drain flows to mass balance
290
2/2
✓ Branch 0 taken 576 times.
✓ Branch 1 taken 3193270 times.
3193846 if ( Subcatch[j].lidArea > 0.0 )
291 {
292 576 lid_addDrainLoads(j, Subcatch[j].newQual, tStep);
293 }
294 }
295
296 //=============================================================================
297
298 3354004 double surfqual_getWtdWashoff(int j, int p, double f)
299 //
300 // Input: j = subcatchment index
301 // p = pollutant index
302 // f = weighting factor
303 // Output: returns pollutant washoff value
304 // Purpose: finds wtd. combination of old and new washoff for a pollutant.
305 //
306 {
307 6708008 return (1.0 - f) * Subcatch[j].oldRunoff * Subcatch[j].oldQual[p] +
308 3354004 f * Subcatch[j].newRunoff *Subcatch[j].newQual[p];
309 }
310
311 //=============================================================================
312
313 3193846 void findPondedLoads(int j, double tStep)
314 //
315 // Input: j = subcatchment index
316 // tStep = time step (sec)
317 // Output: updates pondedQual and OutflowLoad
318 // Purpose: mixes wet deposition and runon pollutant loading with existing
319 // ponded pollutant mass to compute an ouflow loading.
320 //
321 {
322 int p; // pollutant index
323 double cPonded, // ponded concentration (mass/ft3)
324 wPonded, // pollutant mass in ponded water (mass)
325 bmpRemoval, // load reduction by best mgt. practice (mass)
326 vRain, // volume of direct precipitation (ft3)
327 wRain, // wet deposition pollutant load (mass)
328 wRunon, // external runon pollutant load (mass)
329 wInfil, // pollutant load lost to infiltration (mass)
330 wOutflow, // ponded water contribution to runoff load (mass)
331 fullArea, // full subcatchment area (ft2)
332 nonLidArea; // non-LID area (ft2)
333
334 // --- subcatchment and non-LID areas
335
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3193846 times.
3193846 if ( Subcatch[j].area == Subcatch[j].lidArea ) return;
336 3193846 fullArea = Subcatch[j].area;
337 3193846 nonLidArea = fullArea - Subcatch[j].lidArea;
338
339 // --- compute precip. volume over time step (ft3)
340 3193846 vRain = Subcatch[j].rainfall * nonLidArea * tStep;
341
342
2/2
✓ Branch 0 taken 12616450 times.
✓ Branch 1 taken 3193846 times.
15810296 for (p = 0; p < Nobjects[POLLUT]; p++)
343 {
344 // --- update mass balance for wet deposition
345 12616450 wRain = Pollut[p].pptConcen * LperFT3 * vRain;
346 12616450 massbal_updateLoadingTotals(DEPOSITION_LOAD, p, wRain * Pollut[p].mcf);
347
348 // --- surface is dry and has no runon -- add any remaining mass
349 // to overall mass balance's FINAL_LOAD category
350
2/2
✓ Branch 0 taken 3771646 times.
✓ Branch 1 taken 8844804 times.
12616450 if ( Vinflow == 0.0 )
351 {
352 3771646 massbal_updateLoadingTotals(FINAL_LOAD, p,
353 3771646 Subcatch[j].pondedQual[p] * Pollut[p].mcf);
354 3771646 Subcatch[j].pondedQual[p] = 0.0;
355 }
356 else
357 {
358 // --- find concen. of ponded water
359 // (newQual[] temporarily holds runon mass loading)
360 8844804 wRunon = Subcatch[j].newQual[p] * tStep;
361 8844804 wPonded = Subcatch[j].pondedQual[p] + wRain + wRunon;
362 8844804 cPonded = wPonded / Vinflow;
363
364 // --- mass lost to infiltration
365 8844804 wInfil = cPonded * Vinfil;
366
2/2
✓ Branch 0 taken 8844803 times.
✓ Branch 1 taken 1 time.
8844804 wInfil = MIN(wInfil, wPonded);
367 8844804 massbal_updateLoadingTotals(INFIL_LOAD, p, wInfil * Pollut[p].mcf);
368 8844804 wPonded -= wInfil;
369
370 // --- mass lost to runoff
371 8844804 wOutflow = cPonded * Voutflow;
372
1/2
✓ Branch 0 taken 8844804 times.
✗ Branch 1 not taken.
8844804 wOutflow = MIN(wOutflow, wPonded);
373 8844804 wPonded -= wOutflow;
374
375 // --- reduce outflow load by average BMP removal
376 8844804 bmpRemoval = landuse_getAvgBmpEffic(j, p) * wOutflow;
377 8844804 massbal_updateLoadingTotals(BMP_REMOVAL_LOAD, p,
378 8844804 bmpRemoval*Pollut[p].mcf);
379 8844804 wOutflow -= bmpRemoval;
380
381 // --- update ponded mass (using newly computed ponded depth)
382 8844804 Subcatch[j].pondedQual[p] = cPonded * subcatch_getDepth(j) * nonLidArea;
383 8844804 OutflowLoad[p] += wOutflow;
384 }
385 }
386 }
387
388 //=============================================================================
389
390 3193846 void findWashoffLoads(int j, double runoff)
391 //
392 // Input: j = subcatchment index
393 // runoff = subcatchment runoff before internal re-routing or
394 // LID controls (ft/sec)
395 // Output: updates OutflowLoad array
396 // Purpose: computes pollutant washoff loads for each land use and adds these
397 // to the subcatchment's total outflow loads.
398 //
399 {
400 int i, // land use index
401 p, // pollutant index
402 k; // co-pollutant index
403 double w, // co-pollutant load (mass)
404 3193846 area = Subcatch[j].area; // subcatchment area (ft2)
405
406 // --- examine each land use
407
2/2
✓ Branch 0 taken 2720410 times.
✓ Branch 1 taken 473436 times.
3193846 if ( runoff < MIN_RUNOFF ) return;
408
2/2
✓ Branch 0 taken 26163 times.
✓ Branch 1 taken 473436 times.
499599 for (i = 0; i < Nobjects[LANDUSE]; i++)
409 {
410
2/2
✓ Branch 0 taken 16447 times.
✓ Branch 1 taken 9716 times.
26163 if ( Subcatch[j].landFactor[i].fraction > 0.0 )
411 {
412 // --- compute load generated by washoff function
413
2/2
✓ Branch 0 taken 33651 times.
✓ Branch 1 taken 16447 times.
50098 for (p = 0; p < Nobjects[POLLUT]; p++)
414 {
415 33651 OutflowLoad[p] += landuse_getWashoffLoad(
416 33651 i, p, area, Subcatch[j].landFactor, runoff, Voutflow);
417 }
418 }
419 }
420
421 // --- compute contribution from any co-pollutant
422
2/2
✓ Branch 0 taken 1867063 times.
✓ Branch 1 taken 473436 times.
2340499 for (p = 0; p < Nobjects[POLLUT]; p++)
423 {
424 // --- check if pollutant p has a co-pollutant k
425 1867063 k = Pollut[p].coPollut;
426
2/2
✓ Branch 0 taken 12654 times.
✓ Branch 1 taken 1854409 times.
1867063 if ( k >= 0 )
427 {
428 // --- compute addition to washoff from co-pollutant
429 12654 w = Pollut[p].coFraction * OutflowLoad[k];
430
431 // --- add this washoff to buildup mass balance totals
432 // so that things will balance
433 12654 massbal_updateLoadingTotals(BUILDUP_LOAD, p, w * Pollut[p].mcf);
434
435 // --- then also add it to the total washoff load
436 12654 OutflowLoad[p] += w;
437 }
438 }
439 }
440
441 //=============================================================================
442
443 3193846 void findLidLoads(int j, double tStep)
444 //
445 // Input: j = subcatchment index
446 // tStep = time step (sec)
447 // Output: updates OutflowLoad array
448 // Purpose: finds addition to subcatchment pollutant loads from wet deposition
449 // and upstream runon to LID areas.
450 //
451 {
452 int p; // pollutant index
453 int useRunon; // = 1 if LIDs receive upstream runon loads
454 double lidArea, // area occupied by LID units (ft2)
455 vLidRain, // direct precip. falling on LID areas (ft3)
456 wLidRain, // wet deposition pollut. load on LID areas (mass)
457 wLidRunon; // runon pollut. load seen by LID areas (mass)
458
459 // --- find rainfall volume seen by LIDs
460 3193846 lidArea = Subcatch[j].lidArea;
461
2/2
✓ Branch 0 taken 3193270 times.
✓ Branch 1 taken 576 times.
3193846 if ( lidArea == 0.0 ) return;
462 576 vLidRain = Subcatch[j].rainfall * lidArea * tStep;
463
464 // --- use upstream runon load only if LIDs occupy full subcatchment
465 // (for partial LID coverage, runon loads were directed onto non-LID area)
466 576 useRunon = (lidArea == Subcatch[j].area);
467
468
2/2
✓ Branch 0 taken 576 times.
✓ Branch 1 taken 576 times.
1152 for (p = 0; p < Nobjects[POLLUT]; p++)
469 {
470 // --- wet deposition load on LID area
471 576 wLidRain = Pollut[p].pptConcen * vLidRain * LperFT3;
472 576 massbal_updateLoadingTotals(DEPOSITION_LOAD, p, wLidRain * Pollut[p].mcf);
473
474 // --- runon load to LID area from other subcatchments
475
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 576 times.
576 if ( useRunon ) wLidRunon = Subcatch[j].newQual[p] * tStep;
476 576 else wLidRunon = 0.0;
477
478 // --- update total outflow pollutant load (mass)
479 576 OutflowLoad[p] += wLidRain + wLidRunon;
480 }
481 }
482