GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 88.6% 342 / 0 / 386
Functions: 95.5% 21 / 0 / 22
Branches: 70.5% 172 / 0 / 244

subcatch.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // subcatch.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 // Subcatchment runoff functions.
10 //
11 // Update History
12 // ==============
13 // Build 5.1.008:
14 // - Support added for keeping separate track of drain outflows from LIDs.
15 // - Processing of inflow/outflow volumes over a time step was refactored.
16 // - Reported subcatchment runoff includes both surface runoff and LID
17 // drain flows, even though latter can be routed elsewhere.
18 // - Runon now distributed only over non-LID area of a subcatchment, unless
19 // LID covers full area.
20 // - Pollutant buildup and washoff functions were moved to surfqual.c.
21 // Build 5.1.009:
22 // - Runon for full LID subcatchment added to statistical summary.
23 // Build 5.1.010:
24 // - Fixed a bug introduced in 5.1.008 that forgot to include LID
25 // exfiltration as inflow sent to GW routine.
26 // Build 5.1.011:
27 // - Subcatchment percent imperviousness not allowed to exceed 100.
28 // Build 5.1.012:
29 // - Subcatchment bottom elevation used instead of aquifer's when
30 // saving water table value to results file.
31 // Build 5.1.013:
32 // - Rain gage isUsed property now set in subcatch_validate().
33 // - Cumulative impervious and pervious area runoff volumes added
34 // to subcatchment statistics.
35 // - Support added for monthly adjustment of subcatchment's depression
36 // storage, pervious N, and infiltration.
37 // Build 5.1.015:
38 // - Support added for multiple infiltration methods within a project.
39 // - Only pervious area depression storage receives monthly adjustment.
40 //-----------------------------------------------------------------------------
41 #define _CRT_SECURE_NO_DEPRECATE
42
43 #include <math.h>
44 #include <string.h>
45 #include "headers.h"
46 #include "lid.h"
47 #include "odesolve.h"
48
49 //-----------------------------------------------------------------------------
50 // Constants
51 //-----------------------------------------------------------------------------
52 const double MCOEFF = 1.49; // constant in Manning Eq.
53 const double MEXP = 1.6666667; // exponent in Manning Eq.
54 const double ODETOL = 0.0001; // acceptable error for ODE solver
55
56 //-----------------------------------------------------------------------------
57 // Globally shared variables
58 //-----------------------------------------------------------------------------
59 // Volumes (ft3) for a subcatchment over a time step
60 double Vevap; // evaporation
61 double Vpevap; // pervious area evaporation
62 double Vinfil; // non-LID infiltration
63 double Vinflow; // non-LID precip + snowmelt + runon + ponded water
64 double Voutflow; // non-LID runoff to subcatchment's outlet
65 double VlidIn; // impervious area flow to LID units
66 double VlidInfil; // infiltration from LID units
67 double VlidOut; // surface outflow from LID units
68 double VlidDrain; // drain outflow from LID units
69 double VlidReturn; // LID outflow returned to pervious area
70
71 //-----------------------------------------------------------------------------
72 // Locally shared variables
73 //-----------------------------------------------------------------------------
74 static TSubarea* theSubarea; // subarea to which getDdDt() is applied
75 static double Dstore; // monthly adjusted depression storage (ft)
76 static double Alpha; // monthly adjusted runoff coeff.
77 static char *RunoffRoutingWords[] = { w_OUTLET, w_IMPERV, w_PERV, NULL};
78
79 //-----------------------------------------------------------------------------
80 // External functions (declared in funcs.h)
81 //-----------------------------------------------------------------------------
82 // subcatch_readParams (called from parseLine in input.c)
83 // subcatch_readSubareaParams (called from parseLine in input.c)
84 // subcatch_readLanduseParams (called from parseLine in input.c)
85 // subcatch_readInitBuildup (called from parseLine in input.c)
86
87 // subcatch_validate (called from project_validate)
88 // subcatch_initState (called from project_init)
89
90 // subcatch_setOldState (called from runoff_execute)
91 // subcatch_getRunon (called from runoff_execute)
92 // subcatch_addRunon (called from subcatch_getRunon,
93 // lid_addDrainRunon, & runoff_getOutfallRunon)
94 // subcatch_getRunoff (called from runoff_execute)
95 // subcatch_hadRunoff (called from runoff_execute)
96
97 // subcatch_getFracPerv (called from gwater_initState)
98 // subcatch_getStorage (called from massbal_getRunoffError)
99 // subcatch_getDepth (called from findPondedLoads in surfqual.c)
100
101 // subcatch_getWtdOutflow (called from addWetWeatherInflows in routing.c)
102 // subcatch_getResults (called from output_saveSubcatchResults)
103
104 //-----------------------------------------------------------------------------
105 // Function declarations
106 //-----------------------------------------------------------------------------
107 static void getNetPrecip(int j, double* netPrecip, double tStep);
108 static double getSubareaRunoff(int subcatch, int subarea, double area,
109 double rainfall, double evap, double tStep);
110 static double getSubareaInfil(int j, TSubarea* subarea, double precip,
111 double tStep);
112 static double findSubareaRunoff(TSubarea* subarea, double tRunoff);
113 static void updatePondedDepth(TSubarea* subarea, double* tx);
114 static void getDdDt(double t, double* d, double* dddt);
115 static void adjustSubareaParams(int subareaType, int subcatch);
116
117 //=============================================================================
118
119 2393 int subcatch_readParams(int j, char* tok[], int ntoks)
120 //
121 // Input: j = subcatchment index
122 // tok[] = array of string tokens
123 // ntoks = number of tokens
124 // Output: returns an error code
125 // Purpose: reads subcatchment parameters from a tokenized line of input data.
126 //
127 // Data has format:
128 // Name RainGage Outlet Area %Imperv Width Slope CurbLength Snowpack
129 //
130 {
131 int i, k, m;
132 char* id;
133 double x[9];
134
135 // --- check for enough tokens
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2393 times.
2393 if ( ntoks < 8 ) return error_setInpError(ERR_ITEMS, "");
137
138 // --- check that named subcatch exists
139 2393 id = project_findID(SUBCATCH, tok[0]);
140
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2393 times.
2393 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
141
142 // --- check that rain gage exists
143 2393 k = project_findObject(GAGE, tok[1]);
144
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2393 times.
2393 if ( k < 0 ) return error_setInpError(ERR_NAME, tok[1]);
145 2393 x[0] = k;
146
147 // --- check that outlet node or subcatch exists
148 2393 m = project_findObject(NODE, tok[2]);
149 2393 x[1] = m;
150 2393 m = project_findObject(SUBCATCH, tok[2]);
151 2393 x[2] = m;
152
3/4
✓ Branch 0 taken 89 times.
✓ Branch 1 taken 2304 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 89 times.
2393 if ( x[1] < 0.0 && x[2] < 0.0 )
153 return error_setInpError(ERR_NAME, tok[2]);
154
155 // --- read area, %imperv, width, slope, & curb length
156
2/2
✓ Branch 0 taken 11965 times.
✓ Branch 1 taken 2393 times.
14358 for ( i = 3; i < 8; i++)
157 {
158
2/4
✓ Branch 1 taken 11965 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 11965 times.
11965 if ( ! getDouble(tok[i], &x[i]) || x[i] < 0.0 )
159 return error_setInpError(ERR_NUMBER, tok[i]);
160 }
161
162 // --- if snowmelt object named, check that it exists
163 2393 x[8] = -1;
164
2/2
✓ Branch 0 taken 2304 times.
✓ Branch 1 taken 89 times.
2393 if ( ntoks > 8 )
165 {
166 2304 k = project_findObject(SNOWMELT, tok[8]);
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2304 times.
2304 if ( k < 0 ) return error_setInpError(ERR_NAME, tok[8]);
168 2304 x[8] = k;
169 }
170
171 // --- assign input values to subcatch's properties
172 2393 Subcatch[j].ID = id;
173 2393 Subcatch[j].gage = (int)x[0];
174 2393 Subcatch[j].outNode = (int)x[1];
175 2393 Subcatch[j].outSubcatch = (int)x[2];
176 2393 Subcatch[j].area = x[3] / UCF(LANDAREA);
177
2/2
✓ Branch 0 taken 2392 times.
✓ Branch 1 taken 1 time.
2393 Subcatch[j].fracImperv = MIN(x[4], 100.0) / 100.0;
178 2393 Subcatch[j].width = x[5] / UCF(LENGTH);
179 2393 Subcatch[j].slope = x[6] / 100.0;
180 2393 Subcatch[j].curbLength = x[7];
181 2393 Subcatch[j].nPervPattern = -1;
182 2393 Subcatch[j].dStorePattern = -1;
183 2393 Subcatch[j].infilPattern = -1;
184
185 // --- create the snow pack object if it hasn't already been created
186
2/2
✓ Branch 0 taken 2304 times.
✓ Branch 1 taken 89 times.
2393 if ( x[8] >= 0 )
187 {
188
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2304 times.
2304 if ( !snow_createSnowpack(j, (int)x[8]) )
189 return error_setInpError(ERR_MEMORY, "");
190 }
191 2393 return 0;
192 }
193
194 //=============================================================================
195
196 2393 int subcatch_readSubareaParams(char* tok[], int ntoks)
197 //
198 // Input: tok[] = array of string tokens
199 // ntoks = number of tokens
200 // Output: returns an error code
201 // Purpose: reads subcatchment's subarea parameters from a tokenized
202 // line of input data.
203 //
204 // Data has format:
205 // Subcatch Imperv_N Perv_N Imperv_S Perv_S PctZero RouteTo (PctRouted)
206 //
207 {
208 int i, j, k, m;
209 double x[7];
210
211 // --- check for enough tokens
212
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2393 times.
2393 if ( ntoks < 7 ) return error_setInpError(ERR_ITEMS, "");
213
214 // --- check that named subcatch exists
215 2393 j = project_findObject(SUBCATCH, tok[0]);
216
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2393 times.
2393 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
217
218 // --- read in Mannings n, depression storage, & PctZero values
219
2/2
✓ Branch 0 taken 11965 times.
✓ Branch 1 taken 2393 times.
14358 for (i = 0; i < 5; i++)
220 {
221
2/4
✓ Branch 1 taken 11965 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 11965 times.
11965 if ( ! getDouble(tok[i+1], &x[i]) || x[i] < 0.0 )
222 return error_setInpError(ERR_NAME, tok[i+1]);
223 }
224
225 // --- check for valid runoff routing keyword
226 2393 m = findmatch(tok[6], RunoffRoutingWords);
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2393 times.
2393 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[6]);
228
229 // --- get percent routed parameter if present (default is 100)
230 2393 x[5] = m;
231 2393 x[6] = 1.0;
232
2/2
✓ Branch 0 taken 2303 times.
✓ Branch 1 taken 90 times.
2393 if ( ntoks >= 8 )
233 {
234
3/6
✓ Branch 1 taken 2303 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2303 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 2303 times.
2303 if ( ! getDouble(tok[7], &x[6]) || x[6] < 0.0 || x[6] > 100.0 )
235 return error_setInpError(ERR_NUMBER, tok[7]);
236 2303 x[6] /= 100.0;
237 }
238
239 // --- assign input values to each type of subarea
240 2393 Subcatch[j].subArea[IMPERV0].N = x[0];
241 2393 Subcatch[j].subArea[IMPERV1].N = x[0];
242 2393 Subcatch[j].subArea[PERV].N = x[1];
243
244 2393 Subcatch[j].subArea[IMPERV0].dStore = 0.0;
245 2393 Subcatch[j].subArea[IMPERV1].dStore = x[2] / UCF(RAINDEPTH);
246 2393 Subcatch[j].subArea[PERV].dStore = x[3] / UCF(RAINDEPTH);
247
248 2393 Subcatch[j].subArea[IMPERV0].fArea = Subcatch[j].fracImperv * x[4] / 100.0;
249 2393 Subcatch[j].subArea[IMPERV1].fArea = Subcatch[j].fracImperv * (1.0 - x[4] / 100.0);
250 2393 Subcatch[j].subArea[PERV].fArea = (1.0 - Subcatch[j].fracImperv);
251
252 // --- assume that all runoff from each subarea goes to subcatch outlet
253
2/2
✓ Branch 0 taken 7179 times.
✓ Branch 1 taken 2393 times.
9572 for (i = IMPERV0; i <= PERV; i++)
254 {
255 7179 Subcatch[j].subArea[i].routeTo = TO_OUTLET;
256 7179 Subcatch[j].subArea[i].fOutlet = 1.0;
257 }
258
259 // --- modify routing if pervious runoff routed to impervious area
260 // (fOutlet is the fraction of runoff not routed)
261
262 2393 k = (int)x[5];
263
2/2
✓ Branch 0 taken 2379 times.
✓ Branch 1 taken 14 times.
2393 if ( Subcatch[j].fracImperv == 0.0
264
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2371 times.
2393 || Subcatch[j].fracImperv == 1.0 ) k = TO_OUTLET;
265
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2393 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
2393 if ( k == TO_IMPERV && Subcatch[j].fracImperv )
266 {
267 Subcatch[j].subArea[PERV].routeTo = k;
268 Subcatch[j].subArea[PERV].fOutlet = 1.0 - x[6];
269 }
270
271 // --- modify routing if impervious runoff routed to pervious area
272
2/2
✓ Branch 0 taken 2303 times.
✓ Branch 1 taken 90 times.
2393 if ( k == TO_PERV )
273 {
274 2303 Subcatch[j].subArea[IMPERV0].routeTo = k;
275 2303 Subcatch[j].subArea[IMPERV1].routeTo = k;
276 2303 Subcatch[j].subArea[IMPERV0].fOutlet = 1.0 - x[6];
277 2303 Subcatch[j].subArea[IMPERV1].fOutlet = 1.0 - x[6];
278 }
279 2393 return 0;
280 }
281
282 //=============================================================================
283
284 33 int subcatch_readLanduseParams(char* tok[], int ntoks)
285 //
286 // Input: tok[] = array of string tokens
287 // ntoks = number of tokens
288 // Output: returns an error code
289 // Purpose: reads assignment of landuses to subcatchment from a tokenized
290 // line of input data.
291 //
292 // Data has format:
293 // Subcatch landuse percent .... landuse percent
294 //
295 {
296 int j, k, m;
297 double f;
298
299 // --- check for enough tokens
300
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
33 if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
301
302 // --- check that named subcatch exists
303 33 j = project_findObject(SUBCATCH, tok[0]);
304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
33 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
305
306 // --- process each pair of landuse - percent items
307
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 33 times.
66 for ( k = 2; k <= ntoks; k = k+2)
308 {
309 // --- check that named land use exists and is followed by a percent
310 33 m = project_findObject(LANDUSE, tok[k-1]);
311
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
33 if ( m < 0 ) return error_setInpError(ERR_NAME, tok[k-1]);
312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
33 if ( k+1 > ntoks ) return error_setInpError(ERR_ITEMS, "");
313
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 33 times.
33 if ( ! getDouble(tok[k], &f) )
314 return error_setInpError(ERR_NUMBER, tok[k]);
315
316 // --- store land use fraction in subcatch's landFactor property
317 33 Subcatch[j].landFactor[m].fraction = f/100.0;
318 }
319 33 return 0;
320 }
321
322 //=============================================================================
323
324 int subcatch_readInitBuildup(char* tok[], int ntoks)
325 //
326 // Input: tok[] = array of string tokens
327 // ntoks = number of tokens
328 // Output: returns an error code
329 // Purpose: reads initial pollutant buildup on subcatchment from
330 // tokenized line of input data.
331 //
332 // Data has format:
333 // Subcatch pollut initLoad .... pollut initLoad
334 //
335 {
336 int j, k, m;
337 double x;
338
339 // --- check for enough tokens
340 if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
341
342 // --- check that named subcatch exists
343 j = project_findObject(SUBCATCH, tok[0]);
344 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
345
346 // --- process each pair of pollutant - init. load items
347 for ( k = 2; k <= ntoks; k = k+2)
348 {
349 // --- check for valid pollutant name and loading value
350 m = project_findObject(POLLUT, tok[k-1]);
351 if ( m < 0 ) return error_setInpError(ERR_NAME, tok[k-1]);
352 if ( k+1 > ntoks ) return error_setInpError(ERR_ITEMS, "");
353 if ( ! getDouble(tok[k], &x) )
354 return error_setInpError(ERR_NUMBER, tok[k]);
355
356 // --- store loading in subcatch's initBuildup property
357 Subcatch[j].initBuildup[m] = x;
358 }
359 return 0;
360 }
361
362 //=============================================================================
363
364 2393 void subcatch_validate(int j)
365 //
366 // Input: j = subcatchment index
367 // Output: none
368 // Purpose: checks for valid subcatchment input parameters.
369 //
370 {
371 int i;
372 double area;
373 2393 double nonLidArea = Subcatch[j].area;
374
375 // --- check for ambiguous outlet name
376
3/4
✓ Branch 0 taken 2304 times.
✓ Branch 1 taken 89 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2304 times.
2393 if ( Subcatch[j].outNode >= 0 && Subcatch[j].outSubcatch >= 0 )
377 report_writeErrorMsg(ERR_SUBCATCH_OUTLET, Subcatch[j].ID);
378
379 // --- validate subcatchment's groundwater component
380 2393 gwater_validate(j);
381
382 // --- validate placement of LIDs in the subcatchment
383 2393 nonLidArea -= Subcatch[j].lidArea;
384
385 // --- compute alpha (i.e. WCON in old SWMM) for overland flow
386 // NOTE: the area which contributes to alpha for both imperv
387 // subareas w/ and w/o depression storage is the total imperv area.
388
2/2
✓ Branch 0 taken 7179 times.
✓ Branch 1 taken 2393 times.
9572 for (i = IMPERV0; i <= PERV; i++)
389 {
390
2/2
✓ Branch 0 taken 2393 times.
✓ Branch 1 taken 4786 times.
7179 if ( i == PERV )
391 {
392 2393 area = (1.0 - Subcatch[j].fracImperv) * nonLidArea;
393 }
394 else
395 {
396 4786 area = Subcatch[j].fracImperv * nonLidArea;
397 }
398 7179 Subcatch[j].subArea[i].alpha = 0.0;
399
400 //// Possible change to how sub-area width should be assigned. ////
401 //// area = nonLidArea;
402 /////////////////////////////////////////////////////////////////////
403
404
3/4
✓ Branch 0 taken 7137 times.
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 7137 times.
✗ Branch 3 not taken.
7179 if ( area > 0.0 && Subcatch[j].subArea[i].N > 0.0 )
405 {
406 7137 Subcatch[j].subArea[i].alpha = MCOEFF * Subcatch[j].width / area *
407 7137 sqrt(Subcatch[j].slope) / Subcatch[j].subArea[i].N;
408 }
409 }
410
411 // --- set isUsed property of subcatchment's rain gage
412 2393 i = Subcatch[j].gage;
413
1/2
✓ Branch 0 taken 2393 times.
✗ Branch 1 not taken.
2393 if (i >= 0) Gage[i].isUsed = TRUE;
414
415 2393 }
416
417 //=============================================================================
418
419 2393 void subcatch_initState(int j)
420 //
421 // Input: j = subcatchment index
422 // Output: none
423 // Purpose: Initializes the state of a subcatchment.
424 //
425 {
426 int i;
427
428 // --- initialize rainfall, runoff, & snow depth
429 2393 Subcatch[j].rainfall = 0.0;
430 2393 Subcatch[j].oldRunoff = 0.0;
431 2393 Subcatch[j].newRunoff = 0.0;
432 2393 Subcatch[j].oldSnowDepth = 0.0;
433 2393 Subcatch[j].newSnowDepth = 0.0;
434 2393 Subcatch[j].runon = 0.0;
435 2393 Subcatch[j].evapLoss = 0.0;
436 2393 Subcatch[j].infilLoss = 0.0;
437
438 // --- initialize state of infiltration, groundwater, & snow pack objects
439
2/2
✓ Branch 0 taken 2389 times.
✓ Branch 1 taken 4 times.
2393 if ( Subcatch[j].infil == j ) infil_initState(j);
440
2/2
✓ Branch 0 taken 2241 times.
✓ Branch 1 taken 152 times.
2393 if ( Subcatch[j].groundwater ) gwater_initState(j);
441
2/2
✓ Branch 0 taken 2304 times.
✓ Branch 1 taken 89 times.
2393 if ( Subcatch[j].snowpack ) snow_initSnowpack(j);
442
443 // --- initialize state of sub-areas
444
2/2
✓ Branch 0 taken 7179 times.
✓ Branch 1 taken 2393 times.
9572 for (i = IMPERV0; i <= PERV; i++)
445 {
446 7179 Subcatch[j].subArea[i].depth = 0.0;
447 7179 Subcatch[j].subArea[i].inflow = 0.0;
448 7179 Subcatch[j].subArea[i].runoff = 0.0;
449 }
450
451 // --- initialize runoff quality
452 2393 surfqual_initState(j);
453 2393 }
454
455 //=============================================================================
456
457 3309970 void subcatch_setOldState(int j)
458 //
459 // Input: j = subcatchment index
460 // Output: none
461 // Purpose: replaces old state of subcatchment with new state.
462 //
463 {
464 int i;
465 3309970 Subcatch[j].oldRunoff = Subcatch[j].newRunoff;
466 3309970 Subcatch[j].oldSnowDepth = Subcatch[j].newSnowDepth;
467
2/2
✓ Branch 0 taken 9929910 times.
✓ Branch 1 taken 3309970 times.
13239880 for (i = IMPERV0; i <= PERV; i++)
468 {
469 9929910 Subcatch[j].subArea[i].inflow = 0.0;
470 }
471
2/2
✓ Branch 0 taken 12937554 times.
✓ Branch 1 taken 3309970 times.
16247524 for (i = 0; i < Nobjects[POLLUT]; i++)
472 {
473 12937554 Subcatch[j].oldQual[i] = Subcatch[j].newQual[i];
474 12937554 Subcatch[j].newQual[i] = 0.0;
475 }
476 3309970 lid_setOldGroupState(j);
477 3309970 }
478
479 //=============================================================================
480
481 3045254 double subcatch_getFracPerv(int j)
482 //
483 // Purpose: determines what fraction of subcatchment area, including any LID
484 // area, is pervious.
485 // Input: j = subcatchment index
486 // Output: returns fraction of area with pervious cover
487 //
488 {
489 3045254 double fracPerv = 1.0 - Subcatch[j].fracImperv;
490
491
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3045254 times.
3045254 if ( Subcatch[j].lidArea > 0.0 )
492 {
493 fracPerv = (fracPerv * (Subcatch[j].area - Subcatch[j].lidArea) +
494 lid_getPervArea(j)) / Subcatch[j].area;
495 fracPerv = MIN(fracPerv, 1.0);
496 }
497 3045254 return fracPerv;
498 }
499
500 //=============================================================================
501
502 4786 double subcatch_getStorage(int j)
503 //
504 // Input: j = subcatchment index
505 // Output: returns total volume of stored water (ft3)
506 // Purpose: finds total volume of water stored on a subcatchment's surface
507 // and its LIDs at the current time.
508 //
509 {
510 int i;
511 4786 double v = 0.0;
512
513
2/2
✓ Branch 0 taken 14358 times.
✓ Branch 1 taken 4786 times.
19144 for ( i = IMPERV0; i <= PERV; i++)
514 {
515 14358 v += Subcatch[j].subArea[i].depth * Subcatch[j].subArea[i].fArea;
516 }
517 4786 return v * (Subcatch[j].area - Subcatch[j].lidArea) +
518 4786 lid_getStoredVolume(j);
519 }
520
521 //=============================================================================
522
523 3229694 void subcatch_getRunon(int j)
524 //
525 // Input: j = subcatchment index
526 // Output: none
527 // Purpose: Routes runoff from a subcatchment to its outlet subcatchment
528 // or between its subareas.
529 //
530 {
531 int k; // outlet subcatchment index
532 int p; // pollutant index
533 double q; // runon to outlet subcatchment (ft/sec)
534 double q1, q2; // runoff from imperv. areas (ft/sec)
535 double pervArea; // subcatchment pervious area (ft2)
536
537 // --- add previous period's runoff from this subcatchment to the
538 // runon of the outflow subcatchment, if it exists
539 3229694 k = Subcatch[j].outSubcatch;
540 3229694 q = Subcatch[j].oldRunoff;
541
4/4
✓ Branch 0 taken 105019 times.
✓ Branch 1 taken 3124675 times.
✓ Branch 2 taken 882 times.
✓ Branch 3 taken 104137 times.
3229694 if ( k >= 0 && k != j )
542 {
543 882 subcatch_addRunonFlow(k, q);
544
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 882 times.
882 for (p = 0; p < Nobjects[POLLUT]; p++)
545 {
546 Subcatch[k].newQual[p] += q * Subcatch[j].oldQual[p] * LperFT3;
547 }
548 }
549
550 // --- add any LID underdrain flow sent from this subcatchment to
551 // other subcatchments
552
2/2
✓ Branch 0 taken 4443 times.
✓ Branch 1 taken 3225251 times.
3229694 if ( Subcatch[j].lidArea > 0.0 ) lid_addDrainRunon(j);
553
554 // --- add to sub-area inflow any outflow from other subarea in previous period
555 // (NOTE: no transfer of runoff pollutant load, since runoff loads are
556 // based on runoff flow from entire subcatchment.)
557
558 // --- Case 1: imperv --> perv
559
2/2
✓ Branch 0 taken 3227483 times.
✓ Branch 1 taken 2211 times.
3229694 if ( Subcatch[j].fracImperv < 1.0 &&
560
2/2
✓ Branch 0 taken 3114643 times.
✓ Branch 1 taken 112840 times.
3227483 Subcatch[j].subArea[IMPERV0].routeTo == TO_PERV )
561 {
562 // --- add area-wtd. outflow from imperv1 subarea to perv area inflow
563 3114643 q1 = Subcatch[j].subArea[IMPERV0].runoff *
564 3114643 Subcatch[j].subArea[IMPERV0].fArea;
565 3114643 q2 = Subcatch[j].subArea[IMPERV1].runoff *
566 3114643 Subcatch[j].subArea[IMPERV1].fArea;
567 3114643 q = q1 + q2;
568 3114643 Subcatch[j].subArea[PERV].inflow += q *
569 3114643 (1.0 - Subcatch[j].subArea[IMPERV0].fOutlet) /
570 3114643 Subcatch[j].subArea[PERV].fArea;
571 }
572
573 // --- Case 2: perv --> imperv
574
2/2
✓ Branch 0 taken 3228770 times.
✓ Branch 1 taken 924 times.
3229694 if ( Subcatch[j].fracImperv > 0.0 &&
575
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3228770 times.
3228770 Subcatch[j].subArea[PERV].routeTo == TO_IMPERV &&
576 Subcatch[j].subArea[IMPERV1].fArea > 0.0 )
577 {
578 q = Subcatch[j].subArea[PERV].runoff;
579 Subcatch[j].subArea[IMPERV1].inflow +=
580 q * (1.0 - Subcatch[j].subArea[PERV].fOutlet) *
581 Subcatch[j].subArea[PERV].fArea /
582 Subcatch[j].subArea[IMPERV1].fArea;
583 }
584
585 // --- Add any return flow from LID units to pervious subarea
586
4/4
✓ Branch 0 taken 4443 times.
✓ Branch 1 taken 3225251 times.
✓ Branch 2 taken 3795 times.
✓ Branch 3 taken 648 times.
3229694 if ( Subcatch[j].lidArea > 0.0 && Subcatch[j].fracImperv < 1.0 )
587 {
588 3795 pervArea = Subcatch[j].subArea[PERV].fArea *
589 3795 (Subcatch[j].area - Subcatch[j].lidArea);
590 3795 q = lid_getFlowToPerv(j);
591
2/2
✓ Branch 0 taken 3648 times.
✓ Branch 1 taken 147 times.
3795 if ( pervArea > 0.0 )
592 {
593 3648 Subcatch[j].subArea[PERV].inflow += q / pervArea;
594 }
595 }
596 3229694 }
597
598 //=============================================================================
599
600 1268 void subcatch_addRunonFlow(int k, double q)
601 //
602 // Input: k = subcatchment index
603 // q = runon flow rate (cfs) to subcatchment k
604 // Output: none
605 // Purpose: Updates the total runon flow (ft/s) seen by a subcatchment that
606 // receives runon flow from an upstream subcatchment.
607 //
608 {
609 int i;
610 double nonLidArea;
611
612 // --- distribute runoff from upstream subcatchment (in cfs)
613 // uniformly over the non-LID area of current subcatchment (ft/sec)
614
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1268 times.
1268 if ( Subcatch[k].area <= 0.0 ) return;
615 1268 nonLidArea = Subcatch[k].area - Subcatch[k].lidArea;
616
1/2
✓ Branch 0 taken 1268 times.
✗ Branch 1 not taken.
1268 if ( nonLidArea > 0.0 ) q = q / nonLidArea;
617 else q = q / Subcatch[k].area;
618 1268 Subcatch[k].runon += q;
619
620 // --- assign this flow to the 3 types of subareas
621
2/2
✓ Branch 0 taken 3804 times.
✓ Branch 1 taken 1268 times.
5072 for (i = IMPERV0; i <= PERV; i++)
622 {
623 3804 Subcatch[k].subArea[i].inflow += q;
624 }
625 }
626
627 //=============================================================================
628
629 3229694 double subcatch_getRunoff(int j, double tStep)
630 //
631 // Input: j = subcatchment index
632 // tStep = time step (sec)
633 // Output: returns total runoff produced by subcatchment (ft/sec)
634 // Purpose: Computes runoff & new storage depth for subcatchment.
635 //
636 // The 'runoff' value returned by this function is the total runoff
637 // generated (in ft/sec) by the subcatchment before any internal
638 // re-routing is applied. It is used to compute pollutant washoff.
639 //
640 // The 'outflow' value computed here (in cfs) is the surface runoff
641 // that actually leaves the subcatchment after any LID controls are
642 // applied and is saved to Subcatch[j].newRunoff.
643 //
644 {
645 int i; // subarea index
646 double nonLidArea; // non-LID portion of subcatch area (ft2)
647 double area; // sub-area or subcatchment area (ft2)
648 double netPrecip[3]; // subarea net precipitation (ft/sec)
649 double vRain; // rainfall (+ snowfall) volume (ft3)
650 3229694 double vRunon = 0.0; // runon volume from other areas (ft3)
651 3229694 double vOutflow = 0.0; // runoff volume leaving subcatch (ft3)
652 3229694 double runoff = 0.0; // total runoff flow on subcatch (cfs)
653 3229694 double evapRate = 0.0; // potential evaporation rate (ft/sec)
654 double subAreaRunoff; // sub-area runoff rate (cfs)
655 3229694 double vImpervRunoff = 0.0; // impervious area runoff volume (ft3)
656 3229694 double vPervRunoff = 0.0; // pervious area runoff volume (ft3)
657
658 // --- initialize shared water balance variables
659 3229694 Vevap = 0.0;
660 3229694 Vpevap = 0.0;
661 3229694 Vinfil = 0.0;
662 3229694 Voutflow = 0.0;
663 3229694 VlidIn = 0.0;
664 3229694 VlidInfil = 0.0;
665 3229694 VlidOut = 0.0;
666 3229694 VlidDrain = 0.0;
667 3229694 VlidReturn = 0.0;
668
669 // --- find volume of inflow to non-LID portion of subcatchment as existing
670 // ponded water + any runon volume from upstream areas;
671 // rainfall and snowmelt will be added as each sub-area is analyzed
672 3229694 nonLidArea = Subcatch[j].area - Subcatch[j].lidArea;
673 3229694 vRunon = Subcatch[j].runon * tStep * nonLidArea;
674 3229694 Vinflow = vRunon + subcatch_getDepth(j) * nonLidArea;
675
676 // --- find LID runon only if LID occupies full subcatchment
677
2/2
✓ Branch 0 taken 147 times.
✓ Branch 1 taken 3229547 times.
3229694 if ( nonLidArea == 0.0 )
678 147 vRunon = Subcatch[j].runon * tStep * Subcatch[j].area;
679
680 // --- get net precip. (rainfall + snowfall + snowmelt) on the 3 types
681 // of subcatchment sub-areas and update Vinflow with it
682 3229694 getNetPrecip(j, netPrecip, tStep);
683
684 // --- find potential evaporation rate
685
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3229694 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
3229694 if ( Evap.dryOnly && Subcatch[j].rainfall > 0.0 ) evapRate = 0.0;
686 3229694 else evapRate = Evap.rate;
687
688 // --- set monthly infiltration adjustment factor
689 3229694 infil_setInfilFactor(j);
690
691 // --- examine each type of sub-area (impervious w/o depression storage,
692 // impervious w/ depression storage, and pervious)
693
4/4
✓ Branch 0 taken 3229547 times.
✓ Branch 1 taken 147 times.
✓ Branch 2 taken 9688641 times.
✓ Branch 3 taken 3229547 times.
12918335 if ( nonLidArea > 0.0 ) for (i = IMPERV0; i <= PERV; i++)
694 {
695 // --- get runoff from sub-area updating Vevap, Vpevap,
696 // Vinfil & Voutflow)
697 9688641 area = nonLidArea * Subcatch[j].subArea[i].fArea;
698 19377282 Subcatch[j].subArea[i].runoff =
699 9688641 getSubareaRunoff(j, i, area, netPrecip[i], evapRate, tStep);
700 9688641 subAreaRunoff = Subcatch[j].subArea[i].runoff * area;
701
2/2
✓ Branch 0 taken 3229547 times.
✓ Branch 1 taken 6459094 times.
9688641 if (i == PERV) vPervRunoff = subAreaRunoff * tStep;
702 6459094 else vImpervRunoff += subAreaRunoff * tStep;
703 9688641 runoff += subAreaRunoff;
704 }
705
706 // --- evaluate any LID treatment provided (updating Vevap,
707 // Vpevap, VlidInfil, VlidIn, VlidOut, & VlidDrain)
708
2/2
✓ Branch 0 taken 4443 times.
✓ Branch 1 taken 3225251 times.
3229694 if ( Subcatch[j].lidArea > 0.0 )
709 {
710 4443 lid_getRunoff(j, tStep);
711 }
712
713 // --- update groundwater levels & flows if applicable
714
4/4
✓ Branch 0 taken 3137683 times.
✓ Branch 1 taken 92011 times.
✓ Branch 2 taken 3043013 times.
✓ Branch 3 taken 94670 times.
3229694 if ( !IgnoreGwater && Subcatch[j].groundwater )
715 {
716 3043013 gwater_getGroundwater(j, Vpevap, Vinfil+VlidInfil, tStep);
717 }
718
719 // --- save subcatchment's total loss rates (ft/s)
720 3229694 area = Subcatch[j].area;
721 3229694 Subcatch[j].evapLoss = Vevap / tStep / area;
722 3229694 Subcatch[j].infilLoss = (Vinfil + VlidInfil) / tStep / area;
723
724 // --- find net surface runoff volume
725 // (VlidDrain accounts for LID drain flows)
726 3229694 vOutflow = Voutflow // runoff from all non-LID areas
727 3229694 - VlidIn // runoff treated by LID units
728 3229694 + VlidOut; // runoff from LID units
729 3229694 Subcatch[j].newRunoff = vOutflow / tStep;
730
731 // --- obtain external precip. volume (without any snowmelt)
732 3229694 vRain = Subcatch[j].rainfall * tStep * area;
733
734 // --- update the cumulative stats for this subcatchment
735 3229694 stats_updateSubcatchStats(j, vRain, vRunon, Vevap, Vinfil + VlidInfil,
736 vImpervRunoff, vPervRunoff, vOutflow + VlidDrain,
737 3229694 Subcatch[j].newRunoff + VlidDrain/tStep);
738
739 // --- include this subcatchment's contribution to overall flow balance
740 // only if its outlet is a drainage system node
741
4/4
✓ Branch 0 taken 105019 times.
✓ Branch 1 taken 3124675 times.
✓ Branch 2 taken 882 times.
✓ Branch 3 taken 104137 times.
3229694 if ( Subcatch[j].outNode == -1 && Subcatch[j].outSubcatch != j )
742 {
743 882 vOutflow = 0.0;
744 }
745
746 // --- update mass balances
747 3229694 massbal_updateRunoffTotals(RUNOFF_RAINFALL, vRain);
748 3229694 massbal_updateRunoffTotals(RUNOFF_EVAP, Vevap);
749 3229694 massbal_updateRunoffTotals(RUNOFF_INFIL, Vinfil+VlidInfil);
750 3229694 massbal_updateRunoffTotals(RUNOFF_RUNOFF, vOutflow);
751
752 // --- return area-averaged runoff (ft/s)
753 3229694 return runoff / area;
754 }
755
756 //=============================================================================
757
758 3229694 void getNetPrecip(int j, double* netPrecip, double tStep)
759 {
760 //
761 // Purpose: Finds combined rainfall + snowmelt on a subcatchment.
762 // Input: j = subcatchment index
763 // tStep = time step (sec)
764 // Output: netPrecip = rainfall + snowmelt over each type of subarea (ft/s)
765 //
766 int i, k;
767 3229694 double rainfall = 0.0; // rainfall (ft/sec)
768 3229694 double snowfall = 0.0; // snowfall (ft/sec)
769
770 // --- get current rainfall or snowfall from rain gage (in ft/sec)
771 3229694 k = Subcatch[j].gage;
772
1/2
✓ Branch 0 taken 3229694 times.
✗ Branch 1 not taken.
3229694 if ( k >= 0 )
773 {
774 3229694 gage_getPrecip(k, &rainfall, &snowfall);
775 }
776
777 // --- assign total precip. rate to subcatch's rainfall property
778 3229694 Subcatch[j].rainfall = rainfall + snowfall;
779
780 // --- determine net precipitation input (netPrecip) to each sub-area
781
782 // --- if subcatch has a snowpack, then base netPrecip on possible snow melt
783
3/4
✓ Branch 0 taken 3115063 times.
✓ Branch 1 taken 114631 times.
✓ Branch 2 taken 3115063 times.
✗ Branch 3 not taken.
3229694 if ( Subcatch[j].snowpack && !IgnoreSnowmelt )
784 {
785 3115063 Subcatch[j].newSnowDepth =
786 3115063 snow_getSnowMelt(j, rainfall, snowfall, tStep, netPrecip);
787 }
788
789 // --- otherwise netPrecip is just sum of rainfall & snowfall
790 else
791 {
792
2/2
✓ Branch 0 taken 343893 times.
✓ Branch 1 taken 114631 times.
458524 for (i=IMPERV0; i<=PERV; i++) netPrecip[i] = rainfall + snowfall;
793 }
794 3229694 }
795
796 //=============================================================================
797
798 12074498 double subcatch_getDepth(int j)
799 //
800 // Input: j = subcatchment index
801 // Output: returns average depth of ponded water (ft)
802 // Purpose: finds average depth of water over the non-LID portion of a
803 // subcatchment
804 //
805 {
806 int i;
807 double fArea;
808 12074498 double depth = 0.0;
809
810
2/2
✓ Branch 0 taken 36223494 times.
✓ Branch 1 taken 12074498 times.
48297992 for (i = IMPERV0; i <= PERV; i++)
811 {
812 36223494 fArea = Subcatch[j].subArea[i].fArea;
813
2/2
✓ Branch 0 taken 36214488 times.
✓ Branch 1 taken 9006 times.
36223494 if ( fArea > 0.0 ) depth += Subcatch[j].subArea[i].depth * fArea;
814 }
815 12074498 return depth;
816 }
817
818 //=============================================================================
819
820 2083801 double subcatch_getWtdOutflow(int j, double f)
821 //
822 // Input: j = subcatchment index
823 // f = weighting factor.
824 // Output: returns weighted runoff value
825 // Purpose: computes wtd. combination of old and new subcatchment runoff.
826 //
827 {
828
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2083801 times.
2083801 if ( Subcatch[j].area == 0.0 ) return 0.0;
829 2083801 return (1.0 - f) * Subcatch[j].oldRunoff + f * Subcatch[j].newRunoff;
830 }
831
832 //=============================================================================
833
834 232220 void subcatch_getResults(int j, double f, float x[])
835 //
836 // Input: j = subcatchment index
837 // f = weighting factor
838 // Output: x = array of results
839 // Purpose: computes wtd. combination of old and new subcatchment results.
840 //
841 {
842 int p; // pollutant index
843 int k; // rain gage index
844 232220 double f1 = 1.0 - f;
845 double z;
846 double runoff;
847 TGroundwater* gw; // ptr. to groundwater object
848
849 // --- retrieve rainfall for current report period
850 232220 k = Subcatch[j].gage;
851
1/2
✓ Branch 0 taken 232220 times.
✗ Branch 1 not taken.
232220 if ( k >= 0 ) x[SUBCATCH_RAINFALL] = (float)Gage[k].reportRainfall;
852 else x[SUBCATCH_RAINFALL] = 0.0f;
853
854 // --- retrieve snow depth
855 464440 z = ( f1 * Subcatch[j].oldSnowDepth +
856 232220 f * Subcatch[j].newSnowDepth ) * UCF(RAINDEPTH);
857 232220 x[SUBCATCH_SNOWDEPTH] = (float)z;
858
859 // --- retrieve runoff and losses
860 232220 x[SUBCATCH_EVAP] = (float)(Subcatch[j].evapLoss * UCF(EVAPRATE));
861 232220 x[SUBCATCH_INFIL] = (float)(Subcatch[j].infilLoss * UCF(RAINFALL));
862 232220 runoff = f1 * Subcatch[j].oldRunoff + f * Subcatch[j].newRunoff;
863
864 // --- add any LID drain flow to reported runoff
865
2/2
✓ Branch 0 taken 2904 times.
✓ Branch 1 taken 229316 times.
232220 if ( Subcatch[j].lidArea > 0.0 )
866 {
867 2904 runoff += f1 * lid_getDrainFlow(j, PREVIOUS) +
868 2904 f * lid_getDrainFlow(j, CURRENT);
869 }
870
871 // --- if runoff is really small, report it as zero
872
2/2
✓ Branch 0 taken 170864 times.
✓ Branch 1 taken 61356 times.
232220 if ( runoff < MIN_RUNOFF * Subcatch[j].area ) runoff = 0.0;
873 232220 x[SUBCATCH_RUNOFF] = (float)(runoff * UCF(FLOW));
874
875 // --- retrieve groundwater results
876 232220 gw = Subcatch[j].groundwater;
877
2/2
✓ Branch 0 taken 126714 times.
✓ Branch 1 taken 105506 times.
232220 if ( gw )
878 {
879 126714 z = (f1 * gw->oldFlow + f * gw->newFlow) * Subcatch[j].area * UCF(FLOW);
880 126714 x[SUBCATCH_GW_FLOW] = (float)z;
881 126714 z = (gw->bottomElev + gw->lowerDepth) * UCF(LENGTH);
882 126714 x[SUBCATCH_GW_ELEV] = (float)z;
883 126714 z = gw->theta;
884 126714 x[SUBCATCH_SOIL_MOIST] = (float)z;
885 }
886 else
887 {
888 105506 x[SUBCATCH_GW_FLOW] = 0.0f;
889 105506 x[SUBCATCH_GW_ELEV] = 0.0f;
890 105506 x[SUBCATCH_SOIL_MOIST] = 0.0f;
891 }
892
893 // --- retrieve pollutant washoff
894
4/4
✓ Branch 0 taken 229916 times.
✓ Branch 1 taken 2304 times.
✓ Branch 2 taken 682630 times.
✓ Branch 3 taken 229916 times.
914850 if ( !IgnoreQuality ) for (p = 0; p < Nobjects[POLLUT]; p++ )
895 {
896
2/2
✓ Branch 0 taken 484592 times.
✓ Branch 1 taken 198038 times.
682630 if ( runoff == 0.0 ) z = 0.0;
897 198038 else z = f1 * Subcatch[j].oldQual[p] + f * Subcatch[j].newQual[p];
898 682630 x[SUBCATCH_WASHOFF+p] = (float)z;
899 }
900 232220 }
901
902
903 //=============================================================================
904 // SUB-AREA METHODS
905 //=============================================================================
906
907 9688641 double getSubareaRunoff(int j, int i, double area, double precip, double evap,
908 double tStep)
909 //
910 // Purpose: computes runoff & losses from a subarea over the current time step.
911 // Input: j = subcatchment index
912 // i = subarea index
913 // area = sub-area area (ft2)
914 // precip = rainfall + snowmelt over subarea (ft/sec)
915 // evap = evaporation (ft/sec)
916 // tStep = time step (sec)
917 // Output: returns runoff rate from the sub-area (cfs);
918 // updates shared variables Vinflow, Vevap, Vpevap, Vinfil & Voutflow.
919 //
920 {
921 double tRunoff; // time over which runoff occurs (sec)
922 double surfMoisture; // surface water available (ft/sec)
923 double surfEvap; // evap. used for surface water (ft/sec)
924 9688641 double infil = 0.0; // infiltration rate (ft/sec)
925 9688641 double runoff = 0.0; // runoff rate (ft/sec)
926 TSubarea* subarea; // pointer to subarea being analyzed
927
928 // --- no runoff if no area
929
2/2
✓ Branch 0 taken 5550 times.
✓ Branch 1 taken 9683091 times.
9688641 if ( area == 0.0 ) return 0.0;
930
931 // --- assign pointer to current subarea
932 9683091 subarea = &Subcatch[j].subArea[i];
933
934 // --- assume runoff occurs over entire time step
935 9683091 tRunoff = tStep;
936
937 // --- determine evaporation loss rate
938 9683091 surfMoisture = subarea->depth / tStep;
939
2/2
✓ Branch 0 taken 6476993 times.
✓ Branch 1 taken 3206098 times.
9683091 surfEvap = MIN(surfMoisture, evap);
940
941 // --- compute infiltration loss rate
942
2/2
✓ Branch 0 taken 3227336 times.
✓ Branch 1 taken 6455755 times.
9683091 if ( i == PERV ) infil = getSubareaInfil(j, subarea, precip, tStep);
943
944 // --- add precip to other subarea inflows
945 9683091 subarea->inflow += precip;
946 9683091 surfMoisture += subarea->inflow;
947
948 // --- update total inflow, evaporation & infiltration volumes
949 9683091 Vinflow += precip * area * tStep;
950 9683091 Vevap += surfEvap * area * tStep;
951
2/2
✓ Branch 0 taken 3227336 times.
✓ Branch 1 taken 6455755 times.
9683091 if ( i == PERV ) Vpevap += Vevap;
952 9683091 Vinfil += infil * area * tStep;
953
954 // --- assign adjusted runoff coeff. & storage to shared variables
955 9683091 Alpha = subarea->alpha;
956 9683091 Dstore = subarea->dStore;
957 9683091 adjustSubareaParams(i, j);
958
959 // --- if losses exceed available moisture then no ponded water remains
960
2/2
✓ Branch 0 taken 6452078 times.
✓ Branch 1 taken 3231013 times.
9683091 if ( surfEvap + infil >= surfMoisture )
961 {
962 6452078 subarea->depth = 0.0;
963 }
964
965 // --- otherwise reduce inflow by losses and update depth
966 // of ponded water and time over which runoff occurs
967 else
968 {
969 3231013 subarea->inflow -= surfEvap + infil;
970 3231013 updatePondedDepth(subarea, &tRunoff);
971 }
972
973 // --- compute runoff based on updated ponded depth
974 9683091 runoff = findSubareaRunoff(subarea, tRunoff);
975
976 // --- compute runoff volume leaving subcatchment for mass balance purposes
977 // (fOutlet is the fraction of this subarea's runoff that goes to the
978 // subcatchment outlet as opposed to another subarea of the subcatchment)
979 9683091 Voutflow += subarea->fOutlet * runoff * area * tStep;
980 9683091 return runoff;
981 }
982
983 //=============================================================================
984
985 3227336 double getSubareaInfil(int j, TSubarea* subarea, double precip, double tStep)
986 //
987 // Purpose: computes infiltration rate at current time step.
988 // Input: j = subcatchment index
989 // subarea = ptr. to a subarea
990 // precip = rainfall + snowmelt over subarea (ft/sec)
991 // tStep = time step (sec)
992 // Output: returns infiltration rate (ft/s)
993 //
994 {
995 3227336 double infil = 0.0; // actual infiltration rate (ft/sec)
996
997 // --- compute infiltration rate
998 3227336 infil = infil_getInfil(j, tStep, precip,
999 subarea->inflow, subarea->depth);
1000
1001 // --- limit infiltration rate by available void space in unsaturated
1002 // zone of any groundwater aquifer
1003
4/4
✓ Branch 0 taken 3137683 times.
✓ Branch 1 taken 89653 times.
✓ Branch 2 taken 3043013 times.
✓ Branch 3 taken 94670 times.
3227336 if ( !IgnoreGwater && Subcatch[j].groundwater )
1004 {
1005
1/2
✓ Branch 0 taken 3043013 times.
✗ Branch 1 not taken.
3043013 infil = MIN(infil,
1006 Subcatch[j].groundwater->maxInfilVol/tStep);
1007 }
1008 3227336 return infil;
1009 }
1010
1011 //=============================================================================
1012
1013 9683091 double findSubareaRunoff(TSubarea* subarea, double tRunoff)
1014 //
1015 // Purpose: computes runoff (ft/s) from subarea after current time step.
1016 // Input: subarea = ptr. to a subarea
1017 // tRunoff = time step over which runoff occurs (sec)
1018 // Output: returns runoff rate (ft/s)
1019 //
1020 {
1021 9683091 double xDepth = subarea->depth - Dstore;
1022 9683091 double runoff = 0.0;
1023
1024
2/2
✓ Branch 0 taken 1673373 times.
✓ Branch 1 taken 8009718 times.
9683091 if ( xDepth > ZERO )
1025 {
1026 // --- case where nonlinear routing is used
1027
1/2
✓ Branch 0 taken 1673373 times.
✗ Branch 1 not taken.
1673373 if ( subarea->N > 0.0 )
1028 {
1029 1673373 runoff = Alpha * pow(xDepth, MEXP);
1030 }
1031
1032 // --- case where no routing is used (Mannings N = 0)
1033 else
1034 {
1035 runoff = xDepth / tRunoff;
1036 subarea->depth = Dstore;
1037 }
1038 }
1039 else
1040 {
1041 8009718 runoff = 0.0;
1042 }
1043 9683091 return runoff;
1044 }
1045
1046 //=============================================================================
1047
1048 3231013 void updatePondedDepth(TSubarea* subarea, double* dt)
1049 //
1050 // Input: subarea = ptr. to a subarea,
1051 // dt = time step (sec)
1052 // Output: dt = time ponded depth is above depression storage (sec)
1053 // Purpose: computes new ponded depth over subarea after current time step.
1054 //
1055 {
1056 3231013 double ix = subarea->inflow; // excess inflow to subarea (ft/sec)
1057 double dx; // depth above depression storage (ft)
1058 3231013 double tx = *dt; // time over which dx > 0 (sec)
1059
1060 // --- see if not enough inflow to fill depression storage (dStore)
1061
2/2
✓ Branch 0 taken 1557619 times.
✓ Branch 1 taken 1673394 times.
3231013 if ( subarea->depth + ix*tx <= Dstore )
1062 {
1063 1557619 subarea->depth += ix * tx;
1064 }
1065
1066 // --- otherwise use the ODE solver to integrate flow depth
1067 else
1068 {
1069 // --- if depth < Dstore then fill up Dstore & reduce time step
1070 1673394 dx = Dstore - subarea->depth;
1071
3/4
✓ Branch 0 taken 3527 times.
✓ Branch 1 taken 1669867 times.
✓ Branch 2 taken 3527 times.
✗ Branch 3 not taken.
1673394 if ( dx > 0.0 && ix > 0.0 )
1072 {
1073 3527 tx -= dx / ix;
1074 3527 subarea->depth = Dstore;
1075 }
1076
1077 // --- now integrate depth over remaining time step tx
1078
2/4
✓ Branch 0 taken 1673394 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1673394 times.
✗ Branch 3 not taken.
1673394 if ( Alpha > 0.0 && tx > 0.0 )
1079 {
1080 1673394 theSubarea = subarea;
1081 1673394 odesolve_integrate(&(subarea->depth), 1, 0, tx, ODETOL, tx,
1082 getDdDt);
1083 }
1084 else
1085 {
1086 if ( tx < 0.0 ) tx = 0.0;
1087 subarea->depth += ix * tx;
1088 }
1089 }
1090
1091 // --- do not allow ponded depth to go negative
1092
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 3231010 times.
3231013 if ( subarea->depth < 0.0 ) subarea->depth = 0.0;
1093
1094 // --- replace original time step with time ponded depth
1095 // is above depression storage
1096 3231013 *dt = tx;
1097 3231013 }
1098
1099 //=============================================================================
1100
1101 10192566 void getDdDt(double t, double* d, double* dddt)
1102 //
1103 // Input: t = current time (not used)
1104 // d = stored depth (ft)
1105 // Output dddt = derivative of d with respect to time
1106 // Purpose: evaluates derivative of stored depth w.r.t. time
1107 // for the subarea whose runoff is being computed.
1108 //
1109 {
1110 10192566 double ix = theSubarea->inflow;
1111 10192566 double rx = *d - Dstore;
1112
2/2
✓ Branch 0 taken 2057 times.
✓ Branch 1 taken 10190509 times.
10192566 if ( rx < 0.0 )
1113 {
1114 2057 rx = 0.0;
1115 }
1116 else
1117 {
1118 10190509 rx = Alpha * pow(rx, MEXP);
1119 }
1120 10192566 *dddt = ix - rx;
1121 10192566 }
1122
1123 //=============================================================================
1124
1125 9683091 void adjustSubareaParams(int i, int j)
1126 //
1127 // Input: i = type of subarea being analyzed
1128 // j = index of current subcatchment being analyzed
1129 // Output adjusted values of module-level variables Dstore & Alpha
1130 // Purpose: adjusts a pervious subarea's depression storage and its
1131 // runoff coeff. by month of the year.
1132 //
1133 {
1134 int p; // monthly pattern index
1135 int m; // current month of the year
1136 double f; // adjustment factor
1137
1138
2/2
✓ Branch 0 taken 3227336 times.
✓ Branch 1 taken 6455755 times.
9683091 if (i == PERV)
1139 {
1140 // --- depression storage adjustment
1141 3227336 p = Subcatch[j].dStorePattern;
1142
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3227336 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
3227336 if (p >= 0 && Pattern[p].type == MONTHLY_PATTERN)
1143 {
1144 m = datetime_monthOfYear(getDateTime(OldRunoffTime)) - 1;
1145 f = Pattern[p].factor[m];
1146 if (f >= 0.0) Dstore *= f;
1147 }
1148
1149 // --- roughness adjustment to runoff coeff.
1150 3227336 p = Subcatch[j].nPervPattern;
1151
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3227336 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
3227336 if (p >= 0 && Pattern[p].type == MONTHLY_PATTERN)
1152 {
1153 m = datetime_monthOfYear(getDateTime(OldRunoffTime)) - 1;
1154 f = Pattern[p].factor[m];
1155 if (f <= 0.0) Alpha = 0.0;
1156 else Alpha /= f;
1157 }
1158 }
1159 9683091 }
1160