GCC Code Coverage Report

Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 83.6% 495 / 0 / 592
Functions: 100.0% 20 / 0 / 20
Branches: 59.3% 217 / 0 / 366
lidproc.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // lidproc.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 07/13/23 (Build 5.2.4)
7 // Author: L. Rossman
8 //
9 // This module computes the hydrologic performance of an LID (Low Impact
10 // Development) unit at a given point in time.
11 //
12 // Update History
13 // ==============
14 // Build 5.1.007:
15 // - Euler integration now applied to all LID types except Vegetative
16 // Swale which continues to use successive approximation.
17 // - LID layer flux routines were re-written to more accurately model
18 // flooded conditions.
19 // Build 5.1.008:
20 // - MAX_STATE_VARS replaced with MAX_LAYERS.
21 // - Optional soil layer added to Porous Pavement LID.
22 // - Rooftop Disconnection added to types of LIDs.
23 // - Separate accounting of drain flows added.
24 // - Indicator for currently wet LIDs added.
25 // - Detailed reporting procedure fixed.
26 // - Possibile negative head on Bioretention Cell drain avoided.
27 // - Bug in computing flow through Green Roof drainage mat fixed.
28 // Build 5.1.009:
29 // - Fixed typo in net flux rate for vegetative swale LID.
30 // Build 5.1.010:
31 // - New modified version of Green-Ampt used for surface layer infiltration.
32 // Build 5.1.011:
33 // - Re-named STOR_INFIL to STOR_EXFIL and StorageInfil to StorageExfil to
34 // better reflect their meaning.
35 // - Evaporation rates from sub-surface layers reduced by fraction of
36 // surface that is pervious (applies to block paver systems)
37 // - Flux rate routines for LIDs with underdrains modified to produce more
38 // physically meaningful results.
39 // - Reporting of detailed results re-written.
40 // Build 5.1.012:
41 // - Modified upper limit for soil layer percolation.
42 // - Modified upper limit on surface infiltration into rain gardens.
43 // - Modified upper limit on drain flow for LIDs with storage layers.
44 // - Used re-defined wasDry variable for LID reports to fix duplicate lines.
45 // Build 5.1.013:
46 // - Support added for open/closed head levels and multiplier v. head curve
47 // to control underdrain flow.
48 // - Support added for regenerating pavement permeability at fixed intervals.
49 // Build 5.1.014:
50 // - Fixed failure to initialize all LID layer moisture volumes to 0 before
51 // computing LID unit performance in lidproc_getOutflow.
52 // Build 5.2.0:
53 // - Fixed failure to account for effect of Impervious Surface Fraction on
54 // pavement permeability for Permeable Pavement LID
55 // - Fixed units conversion for pavement depth in detailed report file.
56 // Build 5.2.4:
57 // - Modified flux limits in biocellFluxRates, pavementFluxRates and
58 // trenchFluxRates.
59 // - Corrected head calculation in getStorageDrainRate when unit has both
60 // a soil and pavement layer.
61 //-----------------------------------------------------------------------------
62 #define _CRT_SECURE_NO_DEPRECATE
63
64 #include <stdlib.h>
65 #include <string.h>
66 #include <math.h>
67 #include "lid.h"
68 #include "headers.h"
69
70 //-----------------------------------------------------------------------------
71 // Constants
72 //-----------------------------------------------------------------------------
73 #define STOPTOL 0.00328 // integration error tolerance in ft (= 1 mm)
74 #define MINFLOW 2.3e-8 // flow cutoff for dry conditions (= 0.001 in/hr)
75
76 //-----------------------------------------------------------------------------
77 // Enumerations
78 //-----------------------------------------------------------------------------
79 enum LidLayerTypes {
80 SURF, // surface layer
81 SOIL, // soil layer
82 STOR, // storage layer
83 PAVE, // pavement layer
84 DRAIN}; // underdrain system
85
86 enum LidRptVars {
87 SURF_INFLOW, // inflow to surface layer
88 TOTAL_EVAP, // evaporation rate from all layers
89 SURF_INFIL, // infiltration into surface layer
90 PAVE_PERC, // percolation through pavement layer
91 SOIL_PERC, // percolation through soil layer
92 STOR_EXFIL, // exfiltration out of storage layer
93 SURF_OUTFLOW, // outflow from surface layer
94 STOR_DRAIN, // outflow from storage layer
95 SURF_DEPTH, // ponded depth on surface layer
96 PAVE_DEPTH, // water level in pavement layer
97 SOIL_MOIST, // moisture content of soil layer
98 STOR_DEPTH, // water level in storage layer
99 MAX_RPT_VARS};
100
101 //-----------------------------------------------------------------------------
102 // Imported variables
103 //-----------------------------------------------------------------------------
104 extern char HasWetLids; // TRUE if any LIDs are wet (declared in runoff.c)
105
106 //-----------------------------------------------------------------------------
107 // Local Variables
108 //-----------------------------------------------------------------------------
109 static TLidUnit* theLidUnit; // ptr. to a subcatchment's LID unit
110 static TLidProc* theLidProc; // ptr. to a LID process
111
112 static double Tstep; // current time step (sec)
113 static double EvapRate; // evaporation rate (ft/s)
114 static double MaxNativeInfil; // native soil infil. rate limit (ft/s)
115
116 static double SurfaceInflow; // precip. + runon to LID unit (ft/s)
117 static double SurfaceInfil; // infil. rate from surface layer (ft/s)
118 static double SurfaceEvap; // evap. rate from surface layer (ft/s)
119 static double SurfaceOutflow; // outflow from surface layer (ft/s)
120 static double SurfaceVolume; // volume in surface storage (ft)
121
122 static double PaveEvap; // evap. from pavement layer (ft/s)
123 static double PavePerc; // percolation from pavement layer (ft/s)
124 static double PaveVolume; // volume stored in pavement layer (ft)
125
126 static double SoilEvap; // evap. from soil layer (ft/s)
127 static double SoilPerc; // percolation from soil layer (ft/s)
128 static double SoilVolume; // volume in soil/pavement storage (ft)
129
130 static double StorageInflow; // inflow rate to storage layer (ft/s)
131 static double StorageExfil; // exfil. rate from storage layer (ft/s)
132 static double StorageEvap; // evap.rate from storage layer (ft/s)
133 static double StorageDrain; // underdrain flow rate layer (ft/s)
134 static double StorageVolume; // volume in storage layer (ft)
135
136 static double Xold[MAX_LAYERS]; // previous moisture level in LID layers
137
138 //-----------------------------------------------------------------------------
139 // External Functions (declared in lid.h)
140 //-----------------------------------------------------------------------------
141 // lidproc_initWaterBalance (called by lid_initState)
142 // lidproc_getOutflow (called by evalLidUnit in lid.c)
143 // lidproc_saveResults (called by evalLidUnit in lid.c)
144
145 //-----------------------------------------------------------------------------
146 // Local Functions
147 //-----------------------------------------------------------------------------
148 static void barrelFluxRates(double x[], double f[]);
149 static void biocellFluxRates(double x[], double f[]);
150 static void greenRoofFluxRates(double x[], double f[]);
151 static void pavementFluxRates(double x[], double f[]);
152 static void trenchFluxRates(double x[], double f[]);
153 static void swaleFluxRates(double x[], double f[]);
154 static void roofFluxRates(double x[], double f[]);
155
156 static double getSurfaceOutflowRate(double depth);
157 static double getSurfaceOverflowRate(double* surfaceDepth);
158 static double getPavementPermRate(void);
159 static double getSoilPercRate(double theta);
160 static double getStorageExfilRate(void);
161 static double getStorageDrainRate(double storageDepth, double soilTheta,
162 double paveDepth, double surfaceDepth);
163 static double getDrainMatOutflow(double depth);
164 static void getEvapRates(double surfaceVol, double paveVol,
165 double soilVol, double storageVol, double pervFrac);
166
167 static void updateWaterBalance(TLidUnit *lidUnit, double inflow,
168 double evap, double infil, double surfFlow,
169 double drainFlow, double storage);
170
171 static int modpuls_solve(int n, double* x, double* xOld, double* xPrev,
172 double* xMin, double* xMax, double* xTol,
173 double* qOld, double* q, double dt, double omega,
174 void (*derivs)(double*, double*));
175
176 //=============================================================================
177
178 18 void lidproc_initWaterBalance(TLidUnit *lidUnit, double initVol)
179 //
180 // Purpose: initializes the water balance components of a LID unit.
181 // Input: lidUnit = a particular LID unit
182 // initVol = initial water volume stored in the unit (ft)
183 // Output: none
184 //
185 {
186 18 lidUnit->waterBalance.inflow = 0.0;
187 18 lidUnit->waterBalance.evap = 0.0;
188 18 lidUnit->waterBalance.infil = 0.0;
189 18 lidUnit->waterBalance.surfFlow = 0.0;
190 18 lidUnit->waterBalance.drainFlow = 0.0;
191 18 lidUnit->waterBalance.initVol = initVol;
192 18 lidUnit->waterBalance.finalVol = initVol;
193 18 }
194
195 //=============================================================================
196
197 8475 double lidproc_getOutflow(TLidUnit* lidUnit, TLidProc* lidProc, double inflow,
198 double evap, double infil, double maxInfil,
199 double tStep, double* lidEvap,
200 double* lidInfil, double* lidDrain)
201 //
202 // Purpose: computes runoff outflow from a single LID unit.
203 // Input: lidUnit = ptr. to specific LID unit being analyzed
204 // lidProc = ptr. to generic LID process of the LID unit
205 // inflow = runoff rate captured by LID unit (ft/s)
206 // evap = potential evaporation rate (ft/s)
207 // infil = infiltration rate to native soil (ft/s)
208 // maxInfil = max. infiltration rate to native soil (ft/s)
209 // tStep = time step (sec)
210 // Output: lidEvap = evaporation rate for LID unit (ft/s)
211 // lidInfil = infiltration rate for LID unit (ft/s)
212 // lidDrain = drain flow for LID unit (ft/s)
213 // returns surface runoff rate from the LID unit (ft/s)
214 //
215 {
216 int i;
217 double x[MAX_LAYERS]; // layer moisture levels
218 double xOld[MAX_LAYERS]; // work vector
219 double xPrev[MAX_LAYERS]; // work vector
220 double xMin[MAX_LAYERS]; // lower limit on moisture levels
221 double xMax[MAX_LAYERS]; // upper limit on moisture levels
222 double fOld[MAX_LAYERS]; // previously computed flux rates
223 double f[MAX_LAYERS]; // newly computed flux rates
224
225 // convergence tolerance on moisture levels (ft, moisture fraction , ft)
226 8475 double xTol[MAX_LAYERS] = {STOPTOL, STOPTOL, STOPTOL, STOPTOL};
227
228 8475 double omega = 0.0; // integration time weighting
229
230 //... define a pointer to function that computes flux rates through the LID
231 8475 void (*fluxRates) (double *, double *) = NULL;
232
233 //... save references to the LID process and LID unit
234 8475 theLidProc = lidProc;
235 8475 theLidUnit = lidUnit;
236
237 //... save evap, max. infil. & time step to shared variables
238 8475 EvapRate = evap;
239 8475 MaxNativeInfil = maxInfil;
240 8475 Tstep = tStep;
241
242 //... store current moisture levels in vector x
243 8475 x[SURF] = theLidUnit->surfaceDepth;
244 8475 x[SOIL] = theLidUnit->soilMoisture;
245 8475 x[STOR] = theLidUnit->storageDepth;
246 8475 x[PAVE] = theLidUnit->paveDepth;
247
248 //... initialize layer moisture volumes, flux rates and moisture limits
249 8475 SurfaceVolume = 0.0;
250 8475 PaveVolume = 0.0;
251 8475 SoilVolume = 0.0;
252 8475 StorageVolume = 0.0;
253 8475 SurfaceInflow = inflow;
254 8475 SurfaceInfil = 0.0;
255 8475 SurfaceEvap = 0.0;
256 8475 SurfaceOutflow = 0.0;
257 8475 PaveEvap = 0.0;
258 8475 PavePerc = 0.0;
259 8475 SoilEvap = 0.0;
260 8475 SoilPerc = 0.0;
261 8475 StorageInflow = 0.0;
262 8475 StorageExfil = 0.0;
263 8475 StorageEvap = 0.0;
264 8475 StorageDrain = 0.0;
265
2/2
✓ Branch 0 taken 33900 times.
✓ Branch 1 taken 8475 times.
 42375 for (i = 0; i < MAX_LAYERS; i++)
266 {
267 33900 f[i] = 0.0;
268 33900 fOld[i] = theLidUnit->oldFluxRates[i];
269 33900 xMin[i] = 0.0;
270 33900 xMax[i] = BIG;
271 33900 Xold[i] = x[i];
272 }
273
274 //... find Green-Ampt infiltration from surface layer
275
2/2
✓ Branch 0 taken 1152 times.
✓ Branch 1 taken 7323 times.
 8475 if ( theLidProc->lidType == POROUS_PAVEMENT ) SurfaceInfil = 0.0;
276
2/2
✓ Branch 0 taken 3072 times.
✓ Branch 1 taken 4251 times.
 7323 else if ( theLidUnit->soilInfil.Ks > 0.0 )
277 {
278 3072 SurfaceInfil =
279 3072 grnampt_getInfil(&theLidUnit->soilInfil, Tstep,
280 3072 SurfaceInflow, theLidUnit->surfaceDepth,
281 MOD_GREEN_AMPT);
282 }
283 4251 else SurfaceInfil = infil;
284
285 //... set moisture limits for soil & storage layers
286
2/2
✓ Branch 0 taken 4224 times.
✓ Branch 1 taken 4251 times.
 8475 if ( theLidProc->soil.thickness > 0.0 )
287 {
288 4224 xMin[SOIL] = theLidProc->soil.wiltPoint;
289 4224 xMax[SOIL] = theLidProc->soil.porosity;
290 }
291
2/2
✓ Branch 0 taken 1152 times.
✓ Branch 1 taken 7323 times.
 8475 if ( theLidProc->pavement.thickness > 0.0 )
292 {
293 1152 xMax[PAVE] = theLidProc->pavement.thickness;
294 }
295
2/2
✓ Branch 0 taken 5523 times.
✓ Branch 1 taken 2952 times.
 8475 if ( theLidProc->storage.thickness > 0.0 )
296 {
297 5523 xMax[STOR] = theLidProc->storage.thickness;
298 }
299
2/2
✓ Branch 0 taken 1152 times.
✓ Branch 1 taken 7323 times.
 8475 if ( theLidProc->lidType == GREEN_ROOF )
300 {
301 1152 xMax[STOR] = theLidProc->drainMat.thickness;
302 }
303
304 //... determine which flux rate function to use
305
7/8
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 1152 times.
✓ Branch 2 taken 1152 times.
✓ Branch 3 taken 1152 times.
✓ Branch 4 taken 1299 times.
✓ Branch 5 taken 648 times.
✓ Branch 6 taken 1152 times.
✗ Branch 7 not taken.
 8475 switch (theLidProc->lidType)
306 {
307 1920 case BIO_CELL:
308 1920 case RAIN_GARDEN: fluxRates = &biocellFluxRates; break;
309 1152 case GREEN_ROOF: fluxRates = &greenRoofFluxRates; break;
310 1152 case INFIL_TRENCH: fluxRates = &trenchFluxRates; break;
311 1152 case POROUS_PAVEMENT: fluxRates = &pavementFluxRates; break;
312 1299 case RAIN_BARREL: fluxRates = &barrelFluxRates; break;
313 648 case ROOF_DISCON: fluxRates = &roofFluxRates; break;
314 1152 case VEG_SWALE: fluxRates = &swaleFluxRates;
315 1152 omega = 0.5;
316 1152 break;
317 ✗ default: return 0.0;
318 }
319
320 //... update moisture levels and flux rates over the time step
321 8475 i = modpuls_solve(MAX_LAYERS, x, xOld, xPrev, xMin, xMax, xTol,
322 fOld, f, tStep, omega, fluxRates);
323
324 /** For debugging only ********************************************
325 if (i == 0)
326 {
327 fprintf(Frpt.file,
328 "\n WARNING 09: integration failed to converge at %s %s",
329 theDate, theTime);
330 fprintf(Frpt.file,
331 "\n for LID %s placed in subcatchment %s.",
332 theLidProc->ID, theSubcatch->ID);
333 }
334 *******************************************************************/
335
336 //... add any surface overflow to surface outflow
337
3/4
✓ Branch 0 taken 6024 times.
✓ Branch 1 taken 2451 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6024 times.
 8475 if ( theLidProc->surface.canOverflow || theLidUnit->fullWidth == 0.0 )
338 {
339 2451 SurfaceOutflow += getSurfaceOverflowRate(&x[SURF]);
340 }
341
342 //... save updated results
343 8475 theLidUnit->surfaceDepth = x[SURF];
344 8475 theLidUnit->paveDepth = x[PAVE];
345 8475 theLidUnit->soilMoisture = x[SOIL];
346 8475 theLidUnit->storageDepth = x[STOR];
347
2/2
✓ Branch 0 taken 33900 times.
✓ Branch 1 taken 8475 times.
 42375 for (i = 0; i < MAX_LAYERS; i++) theLidUnit->oldFluxRates[i] = f[i];
348
349 //... assign values to LID unit evaporation, infiltration & drain flow
350 8475 *lidEvap = SurfaceEvap + PaveEvap + SoilEvap + StorageEvap;
351 8475 *lidInfil = StorageExfil;
352 8475 *lidDrain = StorageDrain;
353
354 //... return surface outflow (per unit area) from unit
355 8475 return SurfaceOutflow;
356 }
357
358 //=============================================================================
359
360 8475 void lidproc_saveResults(TLidUnit* lidUnit, double ucfRainfall, double ucfRainDepth)
361 //
362 // Purpose: updates the mass balance for an LID unit and saves
363 // current flux rates to the LID report file.
364 // Input: lidUnit = ptr. to LID unit
365 // ucfRainfall = units conversion factor for rainfall rate
366 // ucfDepth = units conversion factor for rainfall depth
367 // Output: none
368 //
369 {
370 double ucf; // units conversion factor
371 double totalEvap; // total evaporation rate (ft/s)
372 double totalVolume; // total volume stored in LID (ft)
373 double rptVars[MAX_RPT_VARS]; // array of reporting variables
374 8475 int isDry = FALSE; // true if current state of LID is dry
375 char timeStamp[TIME_STAMP_SIZE + 1]; // date/time stamp
376 double elapsedHrs; // elapsed hours
377
378 //... find total evap. rate and stored volume
379 8475 totalEvap = SurfaceEvap + PaveEvap + SoilEvap + StorageEvap;
380 8475 totalVolume = SurfaceVolume + PaveVolume + SoilVolume + StorageVolume;
381
382 //... update mass balance totals
383 8475 updateWaterBalance(theLidUnit, SurfaceInflow, totalEvap, StorageExfil,
384 SurfaceOutflow, StorageDrain, totalVolume);
385
386 //... check if dry-weather conditions hold
387
2/2
✓ Branch 0 taken 5807 times.
✓ Branch 1 taken 2668 times.
 8475 if ( SurfaceInflow < MINFLOW &&
388
2/2
✓ Branch 0 taken 5797 times.
✓ Branch 1 taken 10 times.
 5807 SurfaceOutflow < MINFLOW &&
389
2/2
✓ Branch 0 taken 5050 times.
✓ Branch 1 taken 747 times.
 5797 StorageDrain < MINFLOW &&
390
4/4
✓ Branch 0 taken 2249 times.
✓ Branch 1 taken 2801 times.
✓ Branch 2 taken 1429 times.
✓ Branch 3 taken 820 times.
 5050 StorageExfil < MINFLOW &&
391 totalEvap < MINFLOW
392 1429 ) isDry = TRUE;
393
394 //... update status of HasWetLids
395
2/2
✓ Branch 0 taken 7046 times.
✓ Branch 1 taken 1429 times.
 8475 if ( !isDry ) HasWetLids = TRUE;
396
397 //... write results to LID report file
398
2/2
✓ Branch 0 taken 1728 times.
✓ Branch 1 taken 6747 times.
 8475 if ( lidUnit->rptFile )
399 {
400 //... convert rate results to original units (in/hr or mm/hr)
401 1728 ucf = ucfRainfall;
402 1728 rptVars[SURF_INFLOW] = SurfaceInflow*ucf;
403 1728 rptVars[TOTAL_EVAP] = totalEvap*ucf;
404 1728 rptVars[SURF_INFIL] = SurfaceInfil*ucf;
405 1728 rptVars[PAVE_PERC] = PavePerc*ucf;
406 1728 rptVars[SOIL_PERC] = SoilPerc*ucf;
407 1728 rptVars[STOR_EXFIL] = StorageExfil*ucf;
408 1728 rptVars[SURF_OUTFLOW] = SurfaceOutflow*ucf;
409 1728 rptVars[STOR_DRAIN] = StorageDrain*ucf;
410
411 //... convert storage results to original units (in or mm)
412 1728 ucf = ucfRainDepth;
413 1728 rptVars[SURF_DEPTH] = theLidUnit->surfaceDepth*ucf;
414 1728 rptVars[PAVE_DEPTH] = theLidUnit->paveDepth*ucf;
415 1728 rptVars[SOIL_MOIST] = theLidUnit->soilMoisture;
416 1728 rptVars[STOR_DEPTH] = theLidUnit->storageDepth*ucf;
417
418 //... if the current LID state is wet but the previous state was dry
419 // for more than one period then write the saved previous results
420 // to the report file thus marking the end of a dry period
421
2/4
✓ Branch 0 taken 1728 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1728 times.
 1728 if ( !isDry && theLidUnit->rptFile->wasDry > 1)
422 {
423 ✗ fprintf(theLidUnit->rptFile->file, "%s",
424 ✗ theLidUnit->rptFile->results);
425 }
426
427 //... write the current results to a string which is saved between
428 // reporting periods
429 1728 elapsedHrs = NewRunoffTime / 1000.0 / 3600.0;
430 1728 datetime_getTimeStamp(
431 M_D_Y, getDateTime(NewRunoffTime), TIME_STAMP_SIZE, timeStamp);
432 1728 snprintf(theLidUnit->rptFile->results, sizeof(theLidUnit->rptFile->results),
433 "\n%20s\t %8.3f\t %8.3f\t %8.4f\t %8.3f\t %8.3f\t %8.3f\t %8.3f\t"
434 "%8.3f\t %8.3f\t %8.3f\t %8.3f\t %8.3f\t %8.3f",
435 timeStamp, elapsedHrs, rptVars[0], rptVars[1], rptVars[2],
436 rptVars[3], rptVars[4], rptVars[5], rptVars[6], rptVars[7],
437 rptVars[8], rptVars[9], rptVars[10], rptVars[11]);
438
439 //... if the current LID state is dry
440
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1728 times.
 1728 if ( isDry )
441 {
442 //... if the previous state was wet then write the current
443 // results to file marking the start of a dry period
444 ✗ if ( theLidUnit->rptFile->wasDry == 0 )
445 {
446 ✗ fprintf(theLidUnit->rptFile->file, "%s",
447 ✗ theLidUnit->rptFile->results);
448 }
449
450 //... increment the number of successive dry periods
451 ✗ theLidUnit->rptFile->wasDry++;
452 }
453
454 //... if the current LID state is wet
455 else
456 {
457 //... write the current results to the report file
458 1728 fprintf(theLidUnit->rptFile->file, "%s",
459 1728 theLidUnit->rptFile->results);
460
461 //... re-set the number of successive dry periods to 0
462 1728 theLidUnit->rptFile->wasDry = 0;
463 }
464 }
465 8475 }
466
467 //=============================================================================
468
469 648 void roofFluxRates(double x[], double f[])
470 //
471 // Purpose: computes flux rates for roof disconnection.
472 // Input: x = vector of storage levels
473 // Output: f = vector of flux rates
474 //
475 {
476 648 double surfaceDepth = x[SURF];
477
478 648 getEvapRates(surfaceDepth, 0.0, 0.0, 0.0, 1.0);
479 648 SurfaceVolume = surfaceDepth;
480 648 SurfaceInfil = 0.0;
481
2/2
✓ Branch 0 taken 576 times.
✓ Branch 1 taken 72 times.
 648 if ( theLidProc->surface.alpha > 0.0 )
482 576 SurfaceOutflow = getSurfaceOutflowRate(surfaceDepth);
483 72 else getSurfaceOverflowRate(&surfaceDepth);
484
1/2
✓ Branch 1 taken 648 times.
✗ Branch 2 not taken.
 648 StorageDrain = MIN(theLidProc->drain.coeff/UCF(RAINFALL), SurfaceOutflow);
485 648 SurfaceOutflow -= StorageDrain;
486 648 f[SURF] = (SurfaceInflow - SurfaceEvap - StorageDrain - SurfaceOutflow);
487 648 }
488
489 //=============================================================================
490
491 1152 void greenRoofFluxRates(double x[], double f[])
492 //
493 // Purpose: computes flux rates from the layers of a green roof.
494 // Input: x = vector of storage levels
495 // Output: f = vector of flux rates
496 //
497 {
498 // Moisture level variables
499 double surfaceDepth;
500 double soilTheta;
501 double storageDepth;
502
503 // Intermediate variables
504 double availVolume;
505 double maxRate;
506
507 // Green roof properties
508 1152 double soilThickness = theLidProc->soil.thickness;
509 1152 double storageThickness = theLidProc->storage.thickness;
510 1152 double soilPorosity = theLidProc->soil.porosity;
511 1152 double storageVoidFrac = theLidProc->storage.voidFrac;
512 1152 double soilFieldCap = theLidProc->soil.fieldCap;
513 1152 double soilWiltPoint = theLidProc->soil.wiltPoint;
514
515 //... retrieve moisture levels from input vector
516 1152 surfaceDepth = x[SURF];
517 1152 soilTheta = x[SOIL];
518 1152 storageDepth = x[STOR];
519
520 //... convert moisture levels to volumes
521 1152 SurfaceVolume = surfaceDepth * theLidProc->surface.voidFrac;
522 1152 SoilVolume = soilTheta * soilThickness;
523 1152 StorageVolume = storageDepth * storageVoidFrac;
524
525 //... get ET rates
526 1152 availVolume = SoilVolume - soilWiltPoint * soilThickness;
527 1152 getEvapRates(SurfaceVolume, 0.0, availVolume, StorageVolume, 1.0);
528
2/2
✓ Branch 0 taken 224 times.
✓ Branch 1 taken 928 times.
 1152 if ( soilTheta >= soilPorosity ) StorageEvap = 0.0;
529
530 //... soil layer perc rate
531 1152 SoilPerc = getSoilPercRate(soilTheta);
532
533 //... limit perc rate by available water
534 1152 availVolume = (soilTheta - soilFieldCap) * soilThickness;
535
2/2
✓ Branch 0 taken 772 times.
✓ Branch 1 taken 380 times.
 1152 maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;
536
2/2
✓ Branch 0 taken 801 times.
✓ Branch 1 taken 351 times.
 1152 SoilPerc = MIN(SoilPerc, maxRate);
537
2/2
✓ Branch 0 taken 802 times.
✓ Branch 1 taken 350 times.
 1152 SoilPerc = MAX(SoilPerc, 0.0);
538
539 //... storage (drain mat) outflow rate
540 1152 StorageExfil = 0.0;
541 1152 StorageDrain = getDrainMatOutflow(storageDepth);
542
543 //... unit is full
544
3/4
✓ Branch 0 taken 224 times.
✓ Branch 1 taken 928 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 224 times.
 1152 if ( soilTheta >= soilPorosity && storageDepth >= storageThickness )
545 {
546 //... outflow from both layers equals limiting rate
547 ✗ maxRate = MIN(SoilPerc, StorageDrain);
548 ✗ SoilPerc = maxRate;
549 ✗ StorageDrain = maxRate;
550
551 //... adjust inflow rate to soil layer
552 ✗ SurfaceInfil = MIN(SurfaceInfil, maxRate);
553 }
554
555 //... unit not full
556 else
557 {
558 //... limit drainmat outflow by available storage volume
559 1152 maxRate = storageDepth * storageVoidFrac / Tstep - StorageEvap;
560
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
 1152 if ( storageDepth >= storageThickness ) maxRate += SoilPerc;
561
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 maxRate = MAX(maxRate, 0.0);
562
2/2
✓ Branch 0 taken 689 times.
✓ Branch 1 taken 463 times.
 1152 StorageDrain = MIN(StorageDrain, maxRate);
563
564 //... limit soil perc inflow by unused storage volume
565 1152 maxRate = (storageThickness - storageDepth) * storageVoidFrac / Tstep +
566 1152 StorageDrain + StorageEvap;
567
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 SoilPerc = MIN(SoilPerc, maxRate);
568
569 //... adjust surface infil. so soil porosity not exceeded
570 1152 maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
571 1152 SoilPerc + SoilEvap;
572
2/2
✓ Branch 0 taken 928 times.
✓ Branch 1 taken 224 times.
 1152 SurfaceInfil = MIN(SurfaceInfil, maxRate);
573 }
574
575 // ... find surface outflow rate
576 1152 SurfaceOutflow = getSurfaceOutflowRate(surfaceDepth);
577
578 // ... compute overall layer flux rates
579 1152 f[SURF] = (SurfaceInflow - SurfaceEvap - SurfaceInfil - SurfaceOutflow) /
580 1152 theLidProc->surface.voidFrac;
581 1152 f[SOIL] = (SurfaceInfil - SoilEvap - SoilPerc) /
582 1152 theLidProc->soil.thickness;
583 1152 f[STOR] = (SoilPerc - StorageEvap - StorageDrain) /
584 1152 theLidProc->storage.voidFrac;
585 1152 }
586
587 //=============================================================================
588
589 1920 void biocellFluxRates(double x[], double f[])
590 //
591 // Purpose: computes flux rates from the layers of a bio-retention cell LID.
592 // Input: x = vector of storage levels
593 // Output: f = vector of flux rates
594 //
595 {
596 // Moisture level variables
597 double surfaceDepth;
598 double soilTheta;
599 double storageDepth;
600
601 // Intermediate variables
602 double availVolume;
603 double maxRate;
604
605 // LID layer properties
606 1920 double soilThickness = theLidProc->soil.thickness;
607 1920 double soilPorosity = theLidProc->soil.porosity;
608 1920 double soilFieldCap = theLidProc->soil.fieldCap;
609 1920 double soilWiltPoint = theLidProc->soil.wiltPoint;
610 1920 double storageThickness = theLidProc->storage.thickness;
611 1920 double storageVoidFrac = theLidProc->storage.voidFrac;
612
613 //... retrieve moisture levels from input vector
614 1920 surfaceDepth = x[SURF];
615 1920 soilTheta = x[SOIL];
616 1920 storageDepth = x[STOR];
617
618 //... convert moisture levels to volumes
619 1920 SurfaceVolume = surfaceDepth * theLidProc->surface.voidFrac;
620 1920 SoilVolume = soilTheta * soilThickness;
621 1920 StorageVolume = storageDepth * storageVoidFrac;
622
623 //... get ET rates
624 1920 availVolume = SoilVolume - soilWiltPoint * soilThickness;
625 1920 getEvapRates(SurfaceVolume, 0.0, availVolume, StorageVolume, 1.0);
626
2/2
✓ Branch 0 taken 399 times.
✓ Branch 1 taken 1521 times.
 1920 if ( soilTheta >= soilPorosity ) StorageEvap = 0.0;
627
628 //... soil layer perc rate
629 1920 SoilPerc = getSoilPercRate(soilTheta);
630
631 //... limit perc rate by available water
632 1920 availVolume = (soilTheta - soilFieldCap) * soilThickness;
633
2/2
✓ Branch 0 taken 1859 times.
✓ Branch 1 taken 61 times.
 1920 maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;
634
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 SoilPerc = MIN(SoilPerc, maxRate);
635
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 SoilPerc = MAX(SoilPerc, 0.0);
636
637 //... exfiltration rate out of storage layer
638 1920 StorageExfil = getStorageExfilRate();
639
640 //... underdrain flow rate
641 1920 StorageDrain = 0.0;
642
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 if ( theLidProc->drain.coeff > 0.0 )
643 {
644 ✗ StorageDrain = getStorageDrainRate(storageDepth, soilTheta, 0.0,
645 surfaceDepth);
646 }
647
648 //... special case of no storage layer present
649
2/2
✓ Branch 0 taken 1152 times.
✓ Branch 1 taken 768 times.
 1920 if ( storageThickness == 0.0 )
650 {
651 1152 StorageEvap = 0.0;
652
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 maxRate = MIN(SoilPerc, StorageExfil);
653 1152 SoilPerc = maxRate;
654 1152 StorageExfil = maxRate;
655
656 //... limit surface infil. by unused soil volume
657 1152 maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
658 1152 SoilPerc + SoilEvap;
659
2/2
✓ Branch 0 taken 916 times.
✓ Branch 1 taken 236 times.
 1152 SurfaceInfil = MIN(SurfaceInfil, maxRate);
660 }
661
662 else
663 {
664 //... storage & soil layers are full
665
3/4
✓ Branch 0 taken 163 times.
✓ Branch 1 taken 605 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 163 times.
 768 if ( soilTheta >= soilPorosity && storageDepth >= storageThickness )
666 {
667 //... limiting rate is smaller of soil perc and storage outflow
668 ✗ maxRate = StorageExfil + StorageDrain;
669 ✗ if ( SoilPerc < maxRate )
670 {
671 ✗ maxRate = SoilPerc;
672 ✗ if ( maxRate > StorageExfil ) StorageDrain = maxRate - StorageExfil;
673 else
674 {
675 ✗ StorageExfil = maxRate;
676 ✗ StorageDrain = 0.0;
677 }
678 }
679 ✗ else SoilPerc = maxRate;
680
681 //... apply limiting rate to surface infil.
682 ✗ SurfaceInfil = MIN(SurfaceInfil, maxRate);
683 }
684
685 //... either layer not full
686 else
687 {
688 //... limit storage exfiltration by available storage volume
689 768 maxRate = SoilPerc - StorageEvap + storageDepth*storageVoidFrac/Tstep;
690
2/2
✓ Branch 0 taken 163 times.
✓ Branch 1 taken 605 times.
 768 StorageExfil = MIN(StorageExfil, maxRate);
691
1/2
✓ Branch 0 taken 768 times.
✗ Branch 1 not taken.
 768 StorageExfil = MAX(StorageExfil, 0.0);
692
693 //... limit underdrain flow by volume above drain offset
694
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 768 times.
 768 if ( StorageDrain > 0.0 )
695 {
696 ✗ maxRate = -StorageExfil - StorageEvap;
697 ✗ if ( storageDepth >= storageThickness) maxRate += SoilPerc;
698 ✗ if ( theLidProc->drain.offset <= storageDepth )
699 {
700 ✗ maxRate += (storageDepth - theLidProc->drain.offset) *
701 ✗ storageVoidFrac/Tstep;
702 }
703 ✗ maxRate = MAX(maxRate, 0.0);
704 ✗ StorageDrain = MIN(StorageDrain, maxRate);
705 }
706
707 //... limit soil perc by unused storage volume
708 768 maxRate = StorageExfil + StorageDrain + StorageEvap +
709 768 (storageThickness - storageDepth) *
710 768 storageVoidFrac/Tstep;
711
1/2
✓ Branch 0 taken 768 times.
✗ Branch 1 not taken.
 768 SoilPerc = MIN(SoilPerc, maxRate);
712
713 //... limit surface infil. by unused soil volume
714 768 maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
715 768 SoilPerc + SoilEvap;
716
2/2
✓ Branch 0 taken 605 times.
✓ Branch 1 taken 163 times.
 768 SurfaceInfil = MIN(SurfaceInfil, maxRate);
717 }
718 }
719
720 //... find surface layer outflow rate
721 1920 SurfaceOutflow = getSurfaceOutflowRate(surfaceDepth);
722
723 //... compute overall layer flux rates
724 1920 f[SURF] = (SurfaceInflow - SurfaceEvap - SurfaceInfil - SurfaceOutflow) /
725 1920 theLidProc->surface.voidFrac;
726 1920 f[SOIL] = (SurfaceInfil - SoilEvap - SoilPerc) /
727 1920 theLidProc->soil.thickness;
728
2/2
✓ Branch 0 taken 1152 times.
✓ Branch 1 taken 768 times.
 1920 if ( storageThickness == 0.0 ) f[STOR] = 0.0;
729 768 else f[STOR] = (SoilPerc - StorageEvap - StorageExfil - StorageDrain) /
730 768 theLidProc->storage.voidFrac;
731 1920 }
732
733 //=============================================================================
734
735 1152 void trenchFluxRates(double x[], double f[])
736 //
737 // Purpose: computes flux rates from the layers of an infiltration trench LID.
738 // Input: x = vector of storage levels
739 // Output: f = vector of flux rates
740 //
741 {
742 // Moisture level variables
743 double surfaceDepth;
744 double storageDepth;
745
746 // Intermediate variables
747 double availVolume;
748 double maxRate;
749
750 // Storage layer properties
751 1152 double storageThickness = theLidProc->storage.thickness;
752 1152 double storageVoidFrac = theLidProc->storage.voidFrac;
753
754 //... retrieve moisture levels from input vector
755 1152 surfaceDepth = x[SURF];
756 1152 storageDepth = x[STOR];
757
758 //... convert moisture levels to volumes
759 1152 SurfaceVolume = surfaceDepth * theLidProc->surface.voidFrac;
760 1152 SoilVolume = 0.0;
761 1152 StorageVolume = storageDepth * storageVoidFrac;
762
763 //... get ET rates
764 1152 availVolume = (storageThickness - storageDepth) * storageVoidFrac;
765 1152 getEvapRates(SurfaceVolume, 0.0, 0.0, StorageVolume, 1.0);
766
767 //... no storage evap if surface ponded
768
2/2
✓ Branch 0 taken 719 times.
✓ Branch 1 taken 433 times.
 1152 if ( surfaceDepth > 0.0 ) StorageEvap = 0.0;
769
770 //... nominal storage inflow
771 1152 StorageInflow = SurfaceInflow + SurfaceVolume / Tstep;
772
773 //... exfiltration rate out of storage layer
774 1152 StorageExfil = getStorageExfilRate();
775
776 //... underdrain flow rate
777 1152 StorageDrain = 0.0;
778
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
 1152 if ( theLidProc->drain.coeff > 0.0 )
779 {
780 ✗ StorageDrain = getStorageDrainRate(storageDepth, 0.0, 0.0, surfaceDepth);
781 }
782
783 //... limit storage exfiltration by available storage volume
784 1152 maxRate = StorageInflow - StorageEvap + storageDepth*storageVoidFrac/Tstep;
785
2/2
✓ Branch 0 taken 882 times.
✓ Branch 1 taken 270 times.
 1152 StorageExfil = MIN(StorageExfil, maxRate);
786
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 StorageExfil = MAX(StorageExfil, 0.0);
787
788 //... limit underdrain flow by volume above drain offset
789
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
 1152 if ( StorageDrain > 0.0 )
790 {
791 ✗ maxRate = -StorageExfil - StorageEvap;
792 ✗ if (storageDepth >= storageThickness ) maxRate += StorageInflow;
793 ✗ if ( theLidProc->drain.offset <= storageDepth )
794 {
795 ✗ maxRate += (storageDepth - theLidProc->drain.offset) *
796 ✗ storageVoidFrac/Tstep;
797 }
798 ✗ maxRate = MAX(maxRate, 0.0);
799 ✗ StorageDrain = MIN(StorageDrain, maxRate);
800 }
801
802 //... limit storage inflow to not exceed storage layer capacity
803 1152 maxRate = (storageThickness - storageDepth)*storageVoidFrac/Tstep +
804 1152 StorageExfil + StorageEvap + StorageDrain;
805
2/2
✓ Branch 0 taken 432 times.
✓ Branch 1 taken 720 times.
 1152 StorageInflow = MIN(StorageInflow, maxRate);
806
807 //... equate surface infil to storage inflow
808 1152 SurfaceInfil = StorageInflow;
809
810 //... find surface outflow rate
811 1152 SurfaceOutflow = getSurfaceOutflowRate(surfaceDepth);
812
813 // ... find net fluxes for each layer
814 1152 f[SURF] = (SurfaceInflow - SurfaceEvap - StorageInflow - SurfaceOutflow) /
815 1152 theLidProc->surface.voidFrac;;
816 1152 f[STOR] = (StorageInflow - StorageEvap - StorageExfil - StorageDrain) /
817 1152 theLidProc->storage.voidFrac;
818 1152 f[SOIL] = 0.0;
819 1152 }
820
821 //=============================================================================
822
823 1152 void pavementFluxRates(double x[], double f[])
824 //
825 // Purpose: computes flux rates for the layers of a porous pavement LID.
826 // Input: x = vector of storage levels
827 // Output: f = vector of flux rates
828 //
829 {
830 //... Moisture level variables
831 double surfaceDepth;
832 double paveDepth;
833 double soilTheta;
834 double storageDepth;
835
836 //... Intermediate variables
837 1152 double pervFrac = (1.0 - theLidProc->pavement.impervFrac);
838 double storageInflow; // inflow rate to storage layer (ft/s)
839 double availVolume;
840 double maxRate;
841
842 //... LID layer properties
843 1152 double paveVoidFrac = theLidProc->pavement.voidFrac * pervFrac;
844 1152 double paveThickness = theLidProc->pavement.thickness;
845 1152 double soilThickness = theLidProc->soil.thickness;
846 1152 double soilPorosity = theLidProc->soil.porosity;
847 1152 double soilFieldCap = theLidProc->soil.fieldCap;
848 1152 double soilWiltPoint = theLidProc->soil.wiltPoint;
849 1152 double storageThickness = theLidProc->storage.thickness;
850 1152 double storageVoidFrac = theLidProc->storage.voidFrac;
851
852 //... retrieve moisture levels from input vector
853 1152 surfaceDepth = x[SURF];
854 1152 paveDepth = x[PAVE];
855 1152 soilTheta = x[SOIL];
856 1152 storageDepth = x[STOR];
857
858 //... convert moisture levels to volumes
859 1152 SurfaceVolume = surfaceDepth * theLidProc->surface.voidFrac;
860 1152 PaveVolume = paveDepth * paveVoidFrac;
861 1152 SoilVolume = soilTheta * soilThickness;
862 1152 StorageVolume = storageDepth * storageVoidFrac;
863
864 //... get ET rates
865 1152 availVolume = SoilVolume - soilWiltPoint * soilThickness;
866 1152 getEvapRates(SurfaceVolume, PaveVolume, availVolume, StorageVolume,
867 pervFrac);
868
869 //... no storage evap if soil or pavement layer saturated
870
3/4
✓ Branch 0 taken 801 times.
✓ Branch 1 taken 351 times.
✓ Branch 2 taken 801 times.
✗ Branch 3 not taken.
 1152 if ( paveDepth >= paveThickness ||
871
2/2
✓ Branch 0 taken 41 times.
✓ Branch 1 taken 760 times.
 801 ( soilThickness > 0.0 && soilTheta >= soilPorosity )
872 392 ) StorageEvap = 0.0;
873
874 //... find nominal rate of surface infiltration into pavement layer
875 1152 SurfaceInfil = SurfaceInflow + (SurfaceVolume / Tstep);
876
877 //... find perc rate out of pavement layer
878 1152 PavePerc = getPavementPermRate() * pervFrac;
879
880 //... surface infiltration can't exceed pavement permeability
881
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 SurfaceInfil = MIN(SurfaceInfil, PavePerc);
882
883 //... limit pavement perc by available water
884 1152 maxRate = PaveVolume/Tstep + SurfaceInfil - PaveEvap;
885
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 maxRate = MAX(maxRate, 0.0);
886
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
 1152 PavePerc = MIN(PavePerc, maxRate);
887
888 //... find soil layer perc rate
889
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 if ( soilThickness > 0.0 )
890 {
891 1152 SoilPerc = getSoilPercRate(soilTheta);
892 1152 availVolume = (soilTheta - soilFieldCap) * soilThickness;
893
2/2
✓ Branch 0 taken 1118 times.
✓ Branch 1 taken 34 times.
 1152 maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;
894
2/2
✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 32 times.
 1152 SoilPerc = MIN(SoilPerc, maxRate);
895
2/2
✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 32 times.
 1152 SoilPerc = MAX(SoilPerc, 0.0);
896 }
897 ✗ else SoilPerc = PavePerc;
898
899 //... exfiltration rate out of storage layer
900 1152 StorageExfil = getStorageExfilRate();
901
902 //... underdrain flow rate
903 1152 StorageDrain = 0.0;
904
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
 1152 if ( theLidProc->drain.coeff > 0.0 )
905 {
906 ✗ StorageDrain = getStorageDrainRate(storageDepth, soilTheta, paveDepth,
907 surfaceDepth);
908 }
909
910 //... check for adjacent saturated layers
911
912 //... no soil layer, pavement & storage layers are full
913
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1152 if ( soilThickness == 0.0 &&
914 ✗ storageDepth >= storageThickness &&
915 paveDepth >= paveThickness )
916 {
917 //... pavement outflow can't exceed storage outflow
918 ✗ maxRate = StorageEvap + StorageDrain + StorageExfil;
919 ✗ if ( PavePerc > maxRate ) PavePerc = maxRate;
920
921 //... storage outflow can't exceed pavement outflow
922 else
923 {
924 //... use up available exfiltration capacity first
925 ✗ StorageExfil = MIN(StorageExfil, PavePerc);
926 ✗ StorageDrain = PavePerc - StorageExfil;
927 }
928
929 //... set soil perc to pavement perc
930 ✗ SoilPerc = PavePerc;
931
932 //... limit surface infil. by pavement perc
933 ✗ SurfaceInfil = MIN(SurfaceInfil, PavePerc);
934 }
935
936 //... pavement, soil & storage layers are full
937
2/4
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1152 times.
 1152 else if ( soilThickness > 0 &&
938 ✗ storageDepth >= storageThickness &&
939 ✗ soilTheta >= soilPorosity &&
940 paveDepth >= paveThickness )
941 {
942 //... find which layer has limiting flux rate
943 ✗ maxRate = StorageExfil + StorageDrain;
944 ✗ if ( SoilPerc < maxRate) maxRate = SoilPerc;
945 ✗ else maxRate = MIN(maxRate, PavePerc);
946
947 //... use up available storage exfiltration capacity first
948 ✗ if ( maxRate > StorageExfil ) StorageDrain = maxRate - StorageExfil;
949 else
950 {
951 ✗ StorageExfil = maxRate;
952 ✗ StorageDrain = 0.0;
953 }
954 ✗ SoilPerc = maxRate;
955 ✗ PavePerc = maxRate;
956
957 //... limit surface infil. by pavement perc
958 ✗ SurfaceInfil = MIN(SurfaceInfil, PavePerc);
959 }
960
961 //... storage & soil layers are full
962
2/4
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1152 times.
 1152 else if ( soilThickness > 0.0 &&
963 ✗ storageDepth >= storageThickness &&
964 soilTheta >= soilPorosity )
965 {
966 //... soil perc can't exceed storage outflow
967 ✗ maxRate = StorageDrain + StorageExfil;
968 ✗ if ( SoilPerc > maxRate ) SoilPerc = maxRate;
969
970 //... storage outflow can't exceed soil perc
971 else
972 {
973 //... use up available exfiltration capacity first
974 ✗ StorageExfil = MIN(StorageExfil, SoilPerc);
975 ✗ StorageDrain = SoilPerc - StorageExfil;
976 }
977 ✗ PavePerc = MIN(PavePerc, SoilPerc);
978
979 //... limit surface infil. by available pavement volume
980 ✗ availVolume = (paveThickness - paveDepth) * paveVoidFrac;
981 ✗ maxRate = availVolume / Tstep + PavePerc + PaveEvap;
982 ✗ SurfaceInfil = MIN(SurfaceInfil, maxRate);
983 }
984
985 //... soil and pavement layers are full
986
3/4
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 351 times.
✓ Branch 3 taken 801 times.
 1152 else if ( soilThickness > 0.0 &&
987
1/2
✓ Branch 0 taken 351 times.
✗ Branch 1 not taken.
 351 paveDepth >= paveThickness &&
988 soilTheta >= soilPorosity )
989 {
990
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 351 times.
 351 PavePerc = MIN(PavePerc, SoilPerc);
991 351 SoilPerc = PavePerc;
992
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 349 times.
 351 SurfaceInfil = MIN(SurfaceInfil,PavePerc);
993
1/2
✓ Branch 0 taken 351 times.
✗ Branch 1 not taken.
 351 maxRate = MAX(StorageVolume / Tstep + SoilPerc - StorageEvap, 0.0);
994
1/2
✓ Branch 0 taken 351 times.
✗ Branch 1 not taken.
 351 StorageExfil = MIN(StorageExfil, maxRate);
995 }
996
997 //... no adjoining layers are full
998 else
999 {
1000 //... limit storage exfiltration by available storage volume
1001 // (if no soil layer, SoilPerc is same as PavePerc)
1002 801 maxRate = SoilPerc - StorageEvap + StorageVolume / Tstep;
1003
2/2
✓ Branch 0 taken 34 times.
✓ Branch 1 taken 767 times.
 801 maxRate = MAX(0.0, maxRate);
1004
2/2
✓ Branch 0 taken 41 times.
✓ Branch 1 taken 760 times.
 801 StorageExfil = MIN(StorageExfil, maxRate);
1005
1006 //... limit underdrain flow by volume above drain offset
1007
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 801 times.
 801 if ( StorageDrain > 0.0 )
1008 {
1009 ✗ maxRate = -StorageExfil - StorageEvap;
1010 ✗ if (storageDepth >= storageThickness ) maxRate += SoilPerc;
1011 ✗ if ( theLidProc->drain.offset <= storageDepth )
1012 {
1013 ✗ maxRate += (storageDepth - theLidProc->drain.offset) *
1014 ✗ storageVoidFrac/Tstep;
1015 }
1016 ✗ maxRate = MAX(maxRate, 0.0);
1017 ✗ StorageDrain = MIN(StorageDrain, maxRate);
1018 }
1019
1020 //... limit soil & pavement outflow by unused storage volume
1021 801 availVolume = (storageThickness - storageDepth) * storageVoidFrac;
1022 801 maxRate = availVolume/Tstep + StorageEvap + StorageDrain + StorageExfil;
1023
1/2
✓ Branch 0 taken 801 times.
✗ Branch 1 not taken.
 801 maxRate = MAX(maxRate, 0.0);
1024
1/2
✓ Branch 0 taken 801 times.
✗ Branch 1 not taken.
 801 if ( soilThickness > 0.0 )
1025 {
1026
1/2
✓ Branch 0 taken 801 times.
✗ Branch 1 not taken.
 801 SoilPerc = MIN(SoilPerc, maxRate);
1027 801 maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
1028 SoilPerc;
1029 }
1030
2/2
✓ Branch 0 taken 759 times.
✓ Branch 1 taken 42 times.
 801 PavePerc = MIN(PavePerc, maxRate);
1031
1032 //... limit surface infil. by available pavement volume
1033 801 availVolume = (paveThickness - paveDepth) * paveVoidFrac;
1034 801 maxRate = availVolume / Tstep + PavePerc + PaveEvap;
1035
2/2
✓ Branch 0 taken 799 times.
✓ Branch 1 taken 2 times.
 801 SurfaceInfil = MIN(SurfaceInfil, maxRate);
1036 }
1037
1038 //... surface outflow
1039 1152 SurfaceOutflow = getSurfaceOutflowRate(surfaceDepth);
1040
1041 //... compute overall layer flux rates
1042 1152 f[SURF] = SurfaceInflow - SurfaceEvap - SurfaceInfil - SurfaceOutflow;
1043 1152 f[PAVE] = (SurfaceInfil - PaveEvap - PavePerc) / paveVoidFrac;
1044
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 if ( theLidProc->soil.thickness > 0.0)
1045 {
1046 1152 f[SOIL] = (PavePerc - SoilEvap - SoilPerc) / soilThickness;
1047 1152 storageInflow = SoilPerc;
1048 }
1049 else
1050 {
1051 ✗ f[SOIL] = 0.0;
1052 ✗ storageInflow = PavePerc;
1053 ✗ SoilPerc = 0.0;
1054 }
1055 1152 f[STOR] = (storageInflow - StorageEvap - StorageExfil - StorageDrain) /
1056 storageVoidFrac;
1057 1152 }
1058
1059 //=============================================================================
1060
1061 2473 void swaleFluxRates(double x[], double f[])
1062 //
1063 // Purpose: computes flux rates from a vegetative swale LID.
1064 // Input: x = vector of storage levels
1065 // Output: f = vector of flux rates
1066 //
1067 {
1068 double depth; // depth of surface water in swale (ft)
1069 double topWidth; // top width of full swale (ft)
1070 double botWidth; // bottom width of swale (ft)
1071 double length; // length of swale (ft)
1072 double surfInflow; // inflow rate to swale (cfs)
1073 double surfWidth; // top width at current water depth (ft)
1074 double surfArea; // surface area of current water depth (ft2)
1075 double flowArea; // x-section flow area (ft2)
1076 double lidArea; // surface area of full swale (ft2)
1077 double hydRadius; // hydraulic radius for current depth (ft)
1078 double slope; // slope of swale side wall (run/rise)
1079 double volume; // swale volume at current water depth (ft3)
1080 double dVdT; // change in volume w.r.t. time (cfs)
1081 double dStore; // depression storage depth (ft)
1082 double xDepth; // depth above depression storage (ft)
1083
1084 //... retrieve state variable from work vector
1085 2473 depth = x[SURF];
1086
1/2
✓ Branch 0 taken 2473 times.
✗ Branch 1 not taken.
 2473 depth = MIN(depth, theLidProc->surface.thickness);
1087
1088 //... depression storage depth
1089 2473 dStore = 0.0;
1090
1091 //... get swale's bottom width
1092 // (0.5 ft minimum to avoid numerical problems)
1093 2473 slope = theLidProc->surface.sideSlope;
1094 2473 topWidth = theLidUnit->fullWidth;
1095
1/2
✓ Branch 0 taken 2473 times.
✗ Branch 1 not taken.
 2473 topWidth = MAX(topWidth, 0.5);
1096 2473 botWidth = topWidth - 2.0 * slope * theLidProc->surface.thickness;
1097
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2473 times.
 2473 if ( botWidth < 0.5 )
1098 {
1099 ✗ botWidth = 0.5;
1100 ✗ slope = 0.5 * (topWidth - 0.5) / theLidProc->surface.thickness;
1101 }
1102
1103 //... swale's length
1104 2473 lidArea = theLidUnit->area;
1105 2473 length = lidArea / topWidth;
1106
1107 //... top width, surface area and flow area of current ponded depth
1108 2473 surfWidth = botWidth + 2.0 * slope * depth;
1109 2473 surfArea = length * surfWidth;
1110 2473 flowArea = (depth * (botWidth + slope * depth)) *
1111 2473 theLidProc->surface.voidFrac;
1112
1113 //... wet volume and effective depth
1114 2473 volume = length * flowArea;
1115
1116 //... surface inflow into swale (cfs)
1117 2473 surfInflow = SurfaceInflow * lidArea;
1118
1119 //... ET rate in cfs
1120 2473 SurfaceEvap = EvapRate * surfArea;
1121
2/2
✓ Branch 0 taken 1573 times.
✓ Branch 1 taken 900 times.
 2473 SurfaceEvap = MIN(SurfaceEvap, volume/Tstep);
1122
1123 //... infiltration rate to native soil in cfs
1124 2473 StorageExfil = SurfaceInfil * surfArea;
1125
1126 //... no surface outflow if depth below depression storage
1127 2473 xDepth = depth - dStore;
1128
2/2
✓ Branch 0 taken 1212 times.
✓ Branch 1 taken 1261 times.
 2473 if ( xDepth <= ZERO ) SurfaceOutflow = 0.0;
1129
1130 //... otherwise compute a surface outflow
1131 else
1132 {
1133 //... modify flow area to remove depression storage,
1134 1261 flowArea -= (dStore * (botWidth + slope * dStore)) *
1135 1261 theLidProc->surface.voidFrac;
1136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1261 times.
 1261 if ( flowArea < ZERO ) SurfaceOutflow = 0.0;
1137 else
1138 {
1139 //... compute hydraulic radius
1140 1261 botWidth = botWidth + 2.0 * dStore * slope;
1141 1261 hydRadius = botWidth + 2.0 * xDepth * sqrt(1.0 + slope*slope);
1142 1261 hydRadius = flowArea / hydRadius;
1143
1144 //... use Manning Eqn. to find outflow rate in cfs
1145 1261 SurfaceOutflow = theLidProc->surface.alpha * flowArea *
1146 1261 pow(hydRadius, 2./3.);
1147 }
1148 }
1149
1150 //... net flux rate (dV/dt) in cfs
1151 2473 dVdT = surfInflow - SurfaceEvap - StorageExfil - SurfaceOutflow;
1152
1153 //... when full, any net positive inflow becomes spillage
1154
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2473 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2473 if ( depth == theLidProc->surface.thickness && dVdT > 0.0 )
1155 {
1156 ✗ SurfaceOutflow += dVdT;
1157 ✗ dVdT = 0.0;
1158 }
1159
1160 //... convert flux rates to ft/s
1161 2473 SurfaceEvap /= lidArea;
1162 2473 StorageExfil /= lidArea;
1163 2473 SurfaceOutflow /= lidArea;
1164 2473 f[SURF] = dVdT / surfArea;
1165 2473 f[SOIL] = 0.0;
1166 2473 f[STOR] = 0.0;
1167
1168 //... assign values to layer volumes
1169 2473 SurfaceVolume = volume / lidArea;
1170 2473 SoilVolume = 0.0;
1171 2473 StorageVolume = 0.0;
1172 2473 }
1173
1174 //=============================================================================
1175
1176 1299 void barrelFluxRates(double x[], double f[])
1177 //
1178 // Purpose: computes flux rates for a rain barrel LID.
1179 // Input: x = vector of storage levels
1180 // Output: f = vector of flux rates
1181 //
1182 {
1183 1299 double storageDepth = x[STOR];
1184 double head;
1185 double maxValue;
1186
1187 //... assign values to layer volumes
1188 1299 SurfaceVolume = 0.0;
1189 1299 SoilVolume = 0.0;
1190 1299 StorageVolume = storageDepth;
1191
1192 //... initialize flows
1193 1299 SurfaceInfil = 0.0;
1194 1299 SurfaceOutflow = 0.0;
1195 1299 StorageDrain = 0.0;
1196
1197 //... compute outflow if time since last rain exceeds drain delay
1198 // (dryTime is updated in lid.evalLidUnit at each time step)
1199
2/2
✓ Branch 0 taken 147 times.
✓ Branch 1 taken 1152 times.
 1299 if ( theLidProc->drain.delay == 0.0 ||
1200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 147 times.
 147 theLidUnit->dryTime >= theLidProc->drain.delay )
1201 {
1202 1152 head = storageDepth - theLidProc->drain.offset;
1203
2/2
✓ Branch 0 taken 814 times.
✓ Branch 1 taken 338 times.
 1152 if ( head > 0.0 )
1204 {
1205 814 StorageDrain = getStorageDrainRate(storageDepth, 0.0, 0.0, 0.0);
1206 814 maxValue = (head/Tstep);
1207
2/2
✓ Branch 0 taken 811 times.
✓ Branch 1 taken 3 times.
 814 StorageDrain = MIN(StorageDrain, maxValue);
1208 }
1209 }
1210
1211 //... limit inflow to available storage
1212 1299 StorageInflow = SurfaceInflow;
1213 1299 maxValue = (theLidProc->storage.thickness - storageDepth) / Tstep +
1214 StorageDrain;
1215
2/2
✓ Branch 0 taken 1058 times.
✓ Branch 1 taken 241 times.
 1299 StorageInflow = MIN(StorageInflow, maxValue);
1216 1299 SurfaceInfil = StorageInflow;
1217
1218 //... assign values to layer flux rates
1219 1299 f[SURF] = SurfaceInflow - StorageInflow;
1220 1299 f[STOR] = StorageInflow - StorageDrain;
1221 1299 f[SOIL] = 0.0;
1222 1299 }
1223
1224 //=============================================================================
1225
1226 5952 double getSurfaceOutflowRate(double depth)
1227 //
1228 // Purpose: computes outflow rate from a LID's surface layer.
1229 // Input: depth = depth of ponded water on surface layer (ft)
1230 // Output: returns outflow from surface layer (ft/s)
1231 //
1232 // Note: this function should not be applied to swales or rain barrels.
1233 //
1234 {
1235 double delta;
1236 double outflow;
1237
1238 //... no outflow if ponded depth below storage depth
1239 5952 delta = depth - theLidProc->surface.thickness;
1240
2/2
✓ Branch 0 taken 5176 times.
✓ Branch 1 taken 776 times.
 5952 if ( delta < 0.0 ) return 0.0;
1241
1242 //... compute outflow from overland flow Manning equation
1243 776 outflow = theLidProc->surface.alpha * pow(delta, 5.0/3.0) *
1244 776 theLidUnit->fullWidth / theLidUnit->area;
1245
2/2
✓ Branch 0 taken 401 times.
✓ Branch 1 taken 375 times.
 776 outflow = MIN(outflow, delta / Tstep);
1246 776 return outflow;
1247 }
1248
1249 //=============================================================================
1250
1251 1152 double getPavementPermRate()
1252 //
1253 // Purpose: computes reduced permeability of a pavement layer due to
1254 // clogging.
1255 // Input: none
1256 // Output: returns the reduced permeability of the pavement layer (ft/s).
1257 //
1258 {
1259 1152 double permReduction = 0.0;
1260 1152 double clogFactor= theLidProc->pavement.clogFactor;
1261 1152 double regenDays = theLidProc->pavement.regenDays;
1262
1263 // ... find permeability reduction due to clogging
1264
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1152 times.
 1152 if ( clogFactor > 0.0 )
1265 {
1266 // ... see if permeability regeneration has occurred
1267 // (regeneration is assumed to reduce the total
1268 // volumetric loading that the pavement has received)
1269 ✗ if ( regenDays > 0.0 )
1270 {
1271 ✗ if ( OldRunoffTime / 1000.0 / SECperDAY >= theLidUnit->nextRegenDay )
1272 {
1273 // ... reduce total volume treated by degree of regeneration
1274 ✗ theLidUnit->volTreated *=
1275 ✗ (1.0 - theLidProc->pavement.regenDegree);
1276
1277 // ... update next day that regenration occurs
1278 ✗ theLidUnit->nextRegenDay += regenDays;
1279 }
1280 }
1281
1282 // ... find permeabiity reduction factor
1283 ✗ permReduction = theLidUnit->volTreated / clogFactor;
1284 ✗ permReduction = MIN(permReduction, 1.0);
1285 }
1286
1287 // ... return the effective pavement permeability
1288 1152 return theLidProc->pavement.kSat * (1.0 - permReduction);
1289 }
1290
1291 //=============================================================================
1292
1293 4224 double getSoilPercRate(double theta)
1294 //
1295 // Purpose: computes percolation rate of water through a LID's soil layer.
1296 // Input: theta = moisture content (fraction)
1297 // Output: returns percolation rate within soil layer (ft/s)
1298 //
1299 {
1300 double delta; // moisture deficit
1301
1302 // ... no percolation if soil moisture <= field capacity
1303
2/2
✓ Branch 0 taken 476 times.
✓ Branch 1 taken 3748 times.
 4224 if ( theta <= theLidProc->soil.fieldCap ) return 0.0;
1304
1305 // ... perc rate = unsaturated hydraulic conductivity
1306 3748 delta = theLidProc->soil.porosity - theta;
1307 3748 return theLidProc->soil.kSat * exp(-delta * theLidProc->soil.kSlope);
1308
1309 }
1310
1311 //=============================================================================
1312
1313 4224 double getStorageExfilRate()
1314 //
1315 // Purpose: computes exfiltration rate from storage zone into
1316 // native soil beneath a LID.
1317 // Input: depth = depth of water storage zone (ft)
1318 // Output: returns infiltration rate (ft/s)
1319 //
1320 {
1321 4224 double infil = 0.0;
1322 4224 double clogFactor = 0.0;
1323
1324
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4224 times.
 4224 if ( theLidProc->storage.kSat == 0.0 ) return 0.0;
1325
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4224 times.
 4224 if ( MaxNativeInfil == 0.0 ) return 0.0;
1326
1327 //... reduction due to clogging
1328 4224 clogFactor = theLidProc->storage.clogFactor;
1329
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4224 times.
 4224 if ( clogFactor > 0.0 )
1330 {
1331 ✗ clogFactor = theLidUnit->waterBalance.inflow / clogFactor;
1332 ✗ clogFactor = MIN(clogFactor, 1.0);
1333 }
1334
1335 //... infiltration rate = storage Ksat reduced by any clogging
1336 4224 infil = theLidProc->storage.kSat * (1.0 - clogFactor);
1337
1338 //... limit infiltration rate by any groundwater-imposed limit
1339
1/2
✓ Branch 0 taken 4224 times.
✗ Branch 1 not taken.
 4224 return MIN(infil, MaxNativeInfil);
1340 }
1341
1342 //=============================================================================
1343
1344 814 double getStorageDrainRate(double storageDepth, double soilTheta,
1345 double paveDepth, double surfaceDepth)
1346 //
1347 // Purpose: computes underdrain flow rate in a LID's storage layer.
1348 // Input: storageDepth = depth of water in storage layer (ft)
1349 // soilTheta = moisture content of soil layer
1350 // paveDepth = effective depth of water in pavement layer (ft)
1351 // surfaceDepth = depth of ponded water on surface layer (ft)
1352 // Output: returns flow in underdrain (ft/s)
1353 //
1354 // Note: drain eqn. is evaluated in user's units.
1355 // Note: head on drain is water depth in storage layer plus the
1356 // layers above it (soil, pavement, and surface in that order)
1357 // minus the drain outlet offset.
1358 {
1359 814 int curve = theLidProc->drain.qCurve;
1360 814 double head = storageDepth;
1361 814 double outflow = 0.0;
1362 814 double paveThickness = theLidProc->pavement.thickness;
1363 814 double soilThickness = theLidProc->soil.thickness;
1364 814 double soilPorosity = theLidProc->soil.porosity;
1365 814 double soilFieldCap = theLidProc->soil.fieldCap;
1366 814 double storageThickness = theLidProc->storage.thickness;
1367
1368 // --- storage layer is full
1369
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 573 times.
 814 if ( storageDepth >= storageThickness )
1370 {
1371 // --- a soil layer exists
1372
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 241 times.
 241 if ( soilThickness > 0.0 )
1373 {
1374 // --- increase head by fraction of soil layer saturated
1375 ✗ if ( soilTheta > soilFieldCap )
1376 {
1377 ✗ head += (soilTheta - soilFieldCap) /
1378 ✗ (soilPorosity - soilFieldCap) * soilThickness;
1379
1380 // --- soil layer is saturated, increase head by water
1381 // depth in layer above it
1382 ✗ if ( soilTheta >= soilPorosity )
1383 {
1384 ✗ if ( paveThickness > 0.0 )
1385 {
1386 ✗ head += paveDepth;
1387 ✗ if ( paveDepth >= paveThickness ) head += surfaceDepth;
1388 }
1389 ✗ else head += surfaceDepth;
1390 }
1391 }
1392 }
1393
1394 // --- no soil layer so increase head by water level in pavement
1395 // layer and possibly surface layer
1396
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 241 times.
 241 else if ( paveThickness > 0.0 )
1397 {
1398 ✗ head += paveDepth;
1399 ✗ if ( paveDepth >= paveThickness ) head += surfaceDepth;
1400 }
1401 }
1402
1403 // --- no outflow if:
1404 // a) no prior outflow and head below open threshold
1405 // b) prior outflow and head below closed threshold
1406
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 574 times.
 814 if ( theLidUnit->oldDrainFlow == 0.0 &&
1407
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 240 times.
 240 head <= theLidProc->drain.hOpen ) return 0.0;
1408
2/2
✓ Branch 0 taken 574 times.
✓ Branch 1 taken 240 times.
 814 if ( theLidUnit->oldDrainFlow > 0.0 &&
1409
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 574 times.
 574 head <= theLidProc->drain.hClose ) return 0.0;
1410
1411 // --- make head relative to drain offset
1412 814 head -= theLidProc->drain.offset;
1413
1414 // --- compute drain outflow from underdrain flow equation in user units
1415 // (head in inches or mm, flow rate in in/hr or mm/hr)
1416
1/2
✓ Branch 0 taken 814 times.
✗ Branch 1 not taken.
 814 if ( head > ZERO )
1417 {
1418 // --- convert head to user units
1419 814 head *= UCF(RAINDEPTH);
1420
1421 // --- compute drain outflow in user units
1422 814 outflow = theLidProc->drain.coeff *
1423 814 pow(head, theLidProc->drain.expon);
1424
1425 // --- apply user-supplied control curve to outflow
1426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 814 times.
 814 if (curve >= 0) outflow *= table_lookup(&Curve[curve], head);
1427
1428 // --- convert outflow to ft/s
1429 814 outflow /= UCF(RAINFALL);
1430 }
1431 814 return outflow;
1432 }
1433
1434 //=============================================================================
1435
1436 1152 double getDrainMatOutflow(double depth)
1437 //
1438 // Purpose: computes flow rate through a green roof's drainage mat.
1439 // Input: depth = depth of water in drainage mat (ft)
1440 // Output: returns flow in drainage mat (ft/s)
1441 //
1442 {
1443 //... default is to pass all inflow
1444 1152 double result = SoilPerc;
1445
1446 //... otherwise use Manning eqn. if its parameters were supplied
1447
1/2
✓ Branch 0 taken 1152 times.
✗ Branch 1 not taken.
 1152 if ( theLidProc->drainMat.alpha > 0.0 )
1448 {
1449 1152 result = theLidProc->drainMat.alpha * pow(depth, 5.0/3.0) *
1450 1152 theLidUnit->fullWidth / theLidUnit->area *
1451 1152 theLidProc->drainMat.voidFrac;
1452 }
1453 1152 return result;
1454 }
1455
1456 //=============================================================================
1457
1458 6024 void getEvapRates(double surfaceVol, double paveVol, double soilVol,
1459 double storageVol, double pervFrac)
1460 //
1461 // Purpose: computes surface, pavement, soil, and storage evaporation rates.
1462 // Input: surfaceVol = volume/area of ponded water on surface layer (ft)
1463 // paveVol = volume/area of water in pavement pores (ft)
1464 // soilVol = volume/area of water in soil (or pavement) pores (ft)
1465 // storageVol = volume/area of water in storage layer (ft)
1466 // pervFrac = fraction of surface layer that is pervious
1467 // Output: none
1468 //
1469 {
1470 double availEvap;
1471
1472 //... surface evaporation flux
1473 6024 availEvap = EvapRate;
1474
2/2
✓ Branch 0 taken 3565 times.
✓ Branch 1 taken 2459 times.
 6024 SurfaceEvap = MIN(availEvap, surfaceVol/Tstep);
1475
2/2
✓ Branch 0 taken 4426 times.
✓ Branch 1 taken 1598 times.
 6024 SurfaceEvap = MAX(0.0, SurfaceEvap);
1476
2/2
✓ Branch 0 taken 3565 times.
✓ Branch 1 taken 2459 times.
 6024 availEvap = MAX(0.0, (availEvap - SurfaceEvap));
1477 6024 availEvap *= pervFrac;
1478
1479 //... no subsurface evap if water is infiltrating
1480
2/2
✓ Branch 0 taken 3089 times.
✓ Branch 1 taken 2935 times.
 6024 if ( SurfaceInfil > 0.0 )
1481 {
1482 3089 PaveEvap = 0.0;
1483 3089 SoilEvap = 0.0;
1484 3089 StorageEvap = 0.0;
1485 }
1486 else
1487 {
1488 //... pavement evaporation flux
1489
2/2
✓ Branch 0 taken 953 times.
✓ Branch 1 taken 1982 times.
 2935 PaveEvap = MIN(availEvap, paveVol / Tstep);
1490
2/2
✓ Branch 0 taken 953 times.
✓ Branch 1 taken 1982 times.
 2935 availEvap = MAX(0.0, (availEvap - PaveEvap));
1491
1492 //... soil evaporation flux
1493
2/2
✓ Branch 0 taken 2930 times.
✓ Branch 1 taken 5 times.
 2935 SoilEvap = MIN(availEvap, soilVol / Tstep);
1494
2/2
✓ Branch 0 taken 2930 times.
✓ Branch 1 taken 5 times.
 2935 availEvap = MAX(0.0, (availEvap - SoilEvap));
1495
1496 //... storage evaporation flux
1497
2/2
✓ Branch 0 taken 2930 times.
✓ Branch 1 taken 5 times.
 2935 StorageEvap = MIN(availEvap, storageVol / Tstep);
1498 }
1499 6024 }
1500
1501 //=============================================================================
1502
1503 2523 double getSurfaceOverflowRate(double* surfaceDepth)
1504 //
1505 // Purpose: finds surface overflow rate from a LID unit.
1506 // Input: surfaceDepth = depth of water stored in surface layer (ft)
1507 // Output: returns the overflow rate (ft/s)
1508 //
1509 {
1510 2523 double delta = *surfaceDepth - theLidProc->surface.thickness;
1511
2/2
✓ Branch 0 taken 2216 times.
✓ Branch 1 taken 307 times.
 2523 if ( delta <= 0.0 ) return 0.0;
1512 307 *surfaceDepth = theLidProc->surface.thickness;
1513 307 return delta * theLidProc->surface.voidFrac / Tstep;
1514 }
1515
1516 //=============================================================================
1517
1518 8475 void updateWaterBalance(TLidUnit *lidUnit, double inflow, double evap,
1519 double infil, double surfFlow, double drainFlow, double storage)
1520 //
1521 // Purpose: updates components of the water mass balance for a LID unit
1522 // over the current time step.
1523 // Input: lidUnit = a particular LID unit
1524 // inflow = runon + rainfall to the LID unit (ft/s)
1525 // evap = evaporation rate from the unit (ft/s)
1526 // infil = infiltration out the bottom of the unit (ft/s)
1527 // surfFlow = surface runoff from the unit (ft/s)
1528 // drainFlow = underdrain flow from the unit
1529 // storage = volume of water stored in the unit (ft)
1530 // Output: none
1531 //
1532 {
1533 8475 lidUnit->volTreated += inflow * Tstep;
1534 8475 lidUnit->waterBalance.inflow += inflow * Tstep;
1535 8475 lidUnit->waterBalance.evap += evap * Tstep;
1536 8475 lidUnit->waterBalance.infil += infil * Tstep;
1537 8475 lidUnit->waterBalance.surfFlow += surfFlow * Tstep;
1538 8475 lidUnit->waterBalance.drainFlow += drainFlow * Tstep;
1539 8475 lidUnit->waterBalance.finalVol = storage;
1540 8475 }
1541
1542 //=============================================================================
1543
1544 8475 int modpuls_solve(int n, double* x, double* xOld, double* xPrev,
1545 double* xMin, double* xMax, double* xTol,
1546 double* qOld, double* q, double dt, double omega,
1547 void (*derivs)(double*, double*))
1548 //
1549 // Purpose: solves system of equations dx/dt = q(x) for x at end of time step
1550 // dt using a modified Puls method.
1551 // Input: n = number of state variables
1552 // x = vector of state variables
1553 // xOld = state variable values at start of time step
1554 // xPrev = state variable values from previous iteration
1555 // xMin = lower limits on state variables
1556 // xMax = upper limits on state variables
1557 // xTol = convergence tolerances on state variables
1558 // qOld = flux rates at start of time step
1559 // q = flux rates at end of time step
1560 // dt = time step (sec)
1561 // omega = time weighting parameter (use 0 for Euler method
1562 // or 0.5 for modified Puls method)
1563 // derivs = pointer to function that computes flux rates q as a
1564 // function of state variables x
1565 // Output: returns number of steps required for convergence (or 0 if
1566 // process doesn't converge)
1567 //
1568 {
1569 int i;
1570 int canStop;
1571 8475 int steps = 1;
1572 8475 int maxSteps = 20;
1573
1574 //... initialize state variable values
1575
2/2
✓ Branch 0 taken 33900 times.
✓ Branch 1 taken 8475 times.
 42375 for (i=0; i<n; i++)
1576 {
1577 33900 xOld[i] = x[i];
1578 33900 xPrev[i] = x[i];
1579 }
1580
1581 //... repeat until convergence achieved
1582
2/2
✓ Branch 0 taken 9796 times.
✓ Branch 1 taken 73 times.
 9869 while (steps < maxSteps)
1583 {
1584 //... compute flux rates for current state levels
1585 9796 canStop = 1;
1586 9796 derivs(x, q);
1587
1588 //... update state levels based on current flux rates
1589
2/2
✓ Branch 0 taken 39184 times.
✓ Branch 1 taken 9796 times.
 48980 for (i=0; i<n; i++)
1590 {
1591 39184 x[i] = xOld[i] + (omega*qOld[i] + (1.0 - omega)*q[i]) * dt;
1592
1/2
✓ Branch 0 taken 39184 times.
✗ Branch 1 not taken.
 39184 x[i] = MIN(x[i], xMax[i]);
1593
2/2
✓ Branch 0 taken 37963 times.
✓ Branch 1 taken 1221 times.
 39184 x[i] = MAX(x[i], xMin[i]);
1594
1595
2/2
✓ Branch 0 taken 9892 times.
✓ Branch 1 taken 29292 times.
 39184 if ( omega > 0.0 &&
1596
2/2
✓ Branch 0 taken 1394 times.
✓ Branch 1 taken 8498 times.
 9892 fabs(x[i] - xPrev[i]) > xTol[i] ) canStop = 0;
1597 39184 xPrev[i] = x[i];
1598 }
1599
1600 //... return if process converges
1601
2/2
✓ Branch 0 taken 8402 times.
✓ Branch 1 taken 1394 times.
 9796 if (canStop) return steps;
1602 1394 steps++;
1603 }
1604
1605 //... no convergence so return 0
1606 73 return 0;
1607 }
1608