qualrout.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // qualrout.c | ||
| 3 | // | ||
| 4 | // Project: EPA SWMM5 | ||
| 5 | // Version: 5.2 | ||
| 6 | // Date: 08/01/22 (Build 5.2.1) | ||
| 7 | // Author: L. Rossman | ||
| 8 | // | ||
| 9 | // Water quality routing functions. | ||
| 10 | // | ||
| 11 | // Update History | ||
| 12 | // ============== | ||
| 13 | // Build 5.1.008: | ||
| 14 | // - Pollutant mass lost to seepage flow added to mass balance totals. | ||
| 15 | // - Pollutant concen. increased when evaporation occurs. | ||
| 16 | // Build 5.1.009: | ||
| 17 | // - Criterion for dry link/storage node changed to avoid concen. blowup. | ||
| 18 | // Build 5.1.010: | ||
| 19 | // - Entire module re-written to be more compact and easier to follow. | ||
| 20 | // - Neglible depth limit replaced with a negligible volume limit. | ||
| 21 | // Build 5.1.015: | ||
| 22 | // - Fixed mass balance issue for empty storage nodes that flood. | ||
| 23 | // Build 5.2.0: | ||
| 24 | // - Support added for flow capture by inlet structures. | ||
| 25 | // - Definition of a dry node/link modified. | ||
| 26 | // Build 5.2.1: | ||
| 27 | // - Dry non-storage nodes now have quality determined by inflow. | ||
| 28 | // - Wet non-storage nodes with no inflow now have no change in quality. | ||
| 29 | //----------------------------------------------------------------------------- | ||
| 30 | #define _CRT_SECURE_NO_DEPRECATE | ||
| 31 | |||
| 32 | #include <stdio.h> | ||
| 33 | #include <stdlib.h> | ||
| 34 | #include <math.h> | ||
| 35 | #include "headers.h" | ||
| 36 | |||
| 37 | //----------------------------------------------------------------------------- | ||
| 38 | // Constants | ||
| 39 | //----------------------------------------------------------------------------- | ||
| 40 | static const double ZeroVolume = 0.0353147; // 1 liter in ft3 | ||
| 41 | static const double ZeroDepth = 0.003281; // 1 mm in ft | ||
| 42 | |||
| 43 | //----------------------------------------------------------------------------- | ||
| 44 | // External functions (declared in funcs.h) | ||
| 45 | //----------------------------------------------------------------------------- | ||
| 46 | // qualrout_init (called by swmm_start) | ||
| 47 | // qualrout_execute (called by routing_execute) | ||
| 48 | |||
| 49 | //----------------------------------------------------------------------------- | ||
| 50 | // Function declarations | ||
| 51 | //----------------------------------------------------------------------------- | ||
| 52 | static void findLinkMassFlow(int i, double tStep); | ||
| 53 | static void findNodeQual(int j); | ||
| 54 | static void findLinkQual(int i, double tStep); | ||
| 55 | static void findSFLinkQual(int i, double qSeep, double fEvap, double tStep); | ||
| 56 | static void findStorageQual(int j, double tStep); | ||
| 57 | static void updateHRT(int j, double v, double q, double tStep); | ||
| 58 | static double getReactedQual(int p, double c, double v1, double tStep); | ||
| 59 | static double getMixedQual(double c, double v1, double wIn, double qIn, | ||
| 60 | double tStep); | ||
| 61 | //============================================================================= | ||
| 62 | |||
| 63 | 57 | void qualrout_init() | |
| 64 | // | ||
| 65 | // Input: none | ||
| 66 | // Output: none | ||
| 67 | // Purpose: initializes water quality concentrations in all nodes and links. | ||
| 68 | // | ||
| 69 | { | ||
| 70 | int i, p, isWet; | ||
| 71 | double c; | ||
| 72 | |||
| 73 |
2/2✓ Branch 0 taken 9267 times.
✓ Branch 1 taken 57 times.
|
9324 | for (i = 0; i < Nobjects[NODE]; i++) |
| 74 | { | ||
| 75 | 9267 | isWet = ( Node[i].newDepth > ZeroDepth ); | |
| 76 |
2/2✓ Branch 0 taken 22107 times.
✓ Branch 1 taken 9267 times.
|
31374 | for (p = 0; p < Nobjects[POLLUT]; p++) |
| 77 | { | ||
| 78 | 22107 | c = 0.0; | |
| 79 |
2/2✓ Branch 0 taken 913 times.
✓ Branch 1 taken 21194 times.
|
22107 | if ( isWet ) c = Pollut[p].initConcen; |
| 80 | 22107 | Node[i].oldQual[p] = c; | |
| 81 | 22107 | Node[i].newQual[p] = c; | |
| 82 | } | ||
| 83 | } | ||
| 84 | |||
| 85 |
2/2✓ Branch 0 taken 9541 times.
✓ Branch 1 taken 57 times.
|
9598 | for (i = 0; i < Nobjects[LINK]; i++) |
| 86 | { | ||
| 87 | 9541 | isWet = ( Link[i].newDepth > ZeroDepth ); | |
| 88 |
2/2✓ Branch 0 taken 23314 times.
✓ Branch 1 taken 9541 times.
|
32855 | for (p = 0; p < Nobjects[POLLUT]; p++) |
| 89 | { | ||
| 90 | 23314 | c = 0.0; | |
| 91 |
2/2✓ Branch 0 taken 961 times.
✓ Branch 1 taken 22353 times.
|
23314 | if ( isWet ) c = Pollut[p].initConcen; |
| 92 | 23314 | Link[i].oldQual[p] = c; | |
| 93 | 23314 | Link[i].newQual[p] = c; | |
| 94 | } | ||
| 95 | } | ||
| 96 | 57 | } | |
| 97 | |||
| 98 | //============================================================================= | ||
| 99 | |||
| 100 | 201935 | void qualrout_execute(double tStep) | |
| 101 | // | ||
| 102 | // Input: tStep = routing time step (sec) | ||
| 103 | // Output: none | ||
| 104 | // Purpose: routes water quality constituents through the drainage | ||
| 105 | // network over the current time step. | ||
| 106 | // | ||
| 107 | { | ||
| 108 | int i, j; | ||
| 109 | double qIn, vAvg; | ||
| 110 | |||
| 111 | // --- find mass flow each link contributes to its downstream node | ||
| 112 |
2/2✓ Branch 1 taken 2758617 times.
✓ Branch 2 taken 201935 times.
|
2960552 | for ( i = 0; i < Nobjects[LINK]; i++ ) findLinkMassFlow(i, tStep); |
| 113 | |||
| 114 | // --- find new water quality concentration at each node | ||
| 115 |
2/2✓ Branch 0 taken 2902606 times.
✓ Branch 1 taken 201935 times.
|
3104541 | for (j = 0; j < Nobjects[NODE]; j++) |
| 116 | { | ||
| 117 | // --- get node inflow and average volume | ||
| 118 | 2902606 | Node[j].qualInflow = Node[j].inflow; | |
| 119 | 2902606 | qIn = Node[j].qualInflow; | |
| 120 | 2902606 | vAvg = (Node[j].oldVolume + Node[j].newVolume) / 2.0; | |
| 121 | |||
| 122 | // --- save inflow concentrations if treatment applied | ||
| 123 |
2/2✓ Branch 0 taken 88564 times.
✓ Branch 1 taken 2814042 times.
|
2902606 | if ( Node[j].treatment ) |
| 124 | { | ||
| 125 |
2/2✓ Branch 0 taken 12008 times.
✓ Branch 1 taken 76556 times.
|
88564 | if ( qIn < ZERO ) qIn = 0.0; |
| 126 | 88564 | treatmnt_setInflow(qIn, Node[j].newQual); | |
| 127 | } | ||
| 128 | |||
| 129 | // --- find new quality at the node | ||
| 130 |
3/4✓ Branch 0 taken 2878459 times.
✓ Branch 1 taken 24147 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2878459 times.
|
2902606 | if ( Node[j].type == STORAGE || Node[j].oldVolume > ZeroVolume ) |
| 131 | { | ||
| 132 | 24147 | findStorageQual(j, tStep); | |
| 133 | } | ||
| 134 | 2878459 | else findNodeQual(j); | |
| 135 | |||
| 136 | // --- apply treatment to new quality values | ||
| 137 |
2/2✓ Branch 0 taken 88564 times.
✓ Branch 1 taken 2814042 times.
|
2902606 | if ( Node[j].treatment ) treatmnt_treat(j, qIn, vAvg, tStep); |
| 138 | } | ||
| 139 | |||
| 140 | // --- find new water quality in each link | ||
| 141 |
2/2✓ Branch 1 taken 2758617 times.
✓ Branch 2 taken 201935 times.
|
2960552 | for ( i = 0; i < Nobjects[LINK]; i++ ) findLinkQual(i, tStep); |
| 142 | 201935 | } | |
| 143 | |||
| 144 | //============================================================================= | ||
| 145 | |||
| 146 | 6086608 | double getMixedQual(double c, double v1, double wIn, double qIn, double tStep) | |
| 147 | // | ||
| 148 | // Input: c = concentration in reactor at start of time step (mass/ft3) | ||
| 149 | // v1 = volume in reactor at start of time step (ft3) | ||
| 150 | // wIn = mass inflow rate (mass/sec) | ||
| 151 | // qIn = flow inflow rate (cfs) | ||
| 152 | // tStep = time step (sec) | ||
| 153 | // Output: returns pollutant concentration at end of time step (mass/ft3) | ||
| 154 | // Purpose: finds pollutant concentration within a completely mixed reactor. | ||
| 155 | // | ||
| 156 | { | ||
| 157 | double vIn, cIn, cMax; | ||
| 158 | |||
| 159 | // --- if no inflow then reactor concentration is unchanged | ||
| 160 |
2/2✓ Branch 0 taken 992054 times.
✓ Branch 1 taken 5094554 times.
|
6086608 | if ( qIn <= ZERO ) return c; |
| 161 | |||
| 162 | // --- compute concentration of any inflow | ||
| 163 | 5094554 | vIn = qIn * tStep; | |
| 164 | 5094554 | cIn = wIn * tStep / vIn; | |
| 165 | |||
| 166 | // --- mixture concen. can't exceed either original or inflow concen. | ||
| 167 |
2/2✓ Branch 0 taken 3963771 times.
✓ Branch 1 taken 1130783 times.
|
5094554 | cMax = MAX(c, cIn); |
| 168 | |||
| 169 | // --- mix inflow with current reactor contents | ||
| 170 | 5094554 | c = (c*v1 + wIn*tStep) / (v1 + vIn); | |
| 171 |
2/2✓ Branch 0 taken 5092160 times.
✓ Branch 1 taken 2394 times.
|
5094554 | c = MIN(c, cMax); |
| 172 |
1/2✓ Branch 0 taken 5094554 times.
✗ Branch 1 not taken.
|
5094554 | c = MAX(c, 0.0); |
| 173 | 5094554 | return c; | |
| 174 | } | ||
| 175 | |||
| 176 | |||
| 177 | //============================================================================= | ||
| 178 | |||
| 179 | 2758617 | void findLinkMassFlow(int i, double tStep) | |
| 180 | // | ||
| 181 | // Input: i = link index | ||
| 182 | // tStep = time step (sec) | ||
| 183 | // Output: none | ||
| 184 | // Purpose: adds constituent mass flow out of link to the total | ||
| 185 | // accumulation at the link's downstream node. | ||
| 186 | // | ||
| 187 | // Note: Node[].newQual[], the accumulator variable, already contains | ||
| 188 | // contributions from runoff and other external inflows from | ||
| 189 | // calculations made in routing_execute(). | ||
| 190 | { | ||
| 191 | int j, p; | ||
| 192 | double qLink, w; | ||
| 193 | |||
| 194 | // --- find inflow to downstream node | ||
| 195 | 2758617 | qLink = Link[i].newFlow; | |
| 196 | |||
| 197 | // --- identify index of downstream node | ||
| 198 | 2758617 | j = Link[i].node2; | |
| 199 |
2/2✓ Branch 0 taken 9817 times.
✓ Branch 1 taken 2748800 times.
|
2758617 | if ( qLink < 0.0 ) j = Link[i].node1; |
| 200 | |||
| 201 | // --- flow rate into downstream node (adjusted for inlet capture) | ||
| 202 | 2758617 | qLink = fabs(qLink); | |
| 203 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2758617 times.
|
2758617 | if (RouteModel != DW) |
| 204 | ✗ | qLink -= inlet_capturedFlow(i); | |
| 205 | |||
| 206 | // --- examine each pollutant | ||
| 207 |
2/2✓ Branch 0 taken 6249431 times.
✓ Branch 1 taken 2758617 times.
|
9008048 | for (p = 0; p < Nobjects[POLLUT]; p++) |
| 208 | { | ||
| 209 | // --- temporarily accumulate inflow load in Node[j].newQual | ||
| 210 | 6249431 | w = qLink * Link[i].oldQual[p]; | |
| 211 | 6249431 | Node[j].newQual[p] += w; | |
| 212 | |||
| 213 | // --- update total load transported by link | ||
| 214 | 6249431 | Link[i].totalLoad[p] += w * tStep; | |
| 215 | } | ||
| 216 | 2758617 | Node[j].qualInflow += qLink; | |
| 217 | 2758617 | } | |
| 218 | |||
| 219 | //============================================================================= | ||
| 220 | |||
| 221 | 2878459 | void findNodeQual(int j) | |
| 222 | // | ||
| 223 | // Input: j = node index | ||
| 224 | // Output: none | ||
| 225 | // Purpose: finds new quality in a node with no storage volume. | ||
| 226 | // | ||
| 227 | { | ||
| 228 | int p; | ||
| 229 | double qNode; | ||
| 230 | |||
| 231 | // --- if there is flow into node then concen. = mass inflow/node flow | ||
| 232 | 2878459 | qNode = Node[j].qualInflow; | |
| 233 |
2/2✓ Branch 0 taken 2529949 times.
✓ Branch 1 taken 348510 times.
|
2878459 | if ( qNode > ZERO ) |
| 234 | { | ||
| 235 |
2/2✓ Branch 0 taken 5596142 times.
✓ Branch 1 taken 2529949 times.
|
8126091 | for (p = 0; p < Nobjects[POLLUT]; p++) |
| 236 | { | ||
| 237 | 5596142 | Node[j].newQual[p] /= qNode; | |
| 238 | } | ||
| 239 | } | ||
| 240 | |||
| 241 | // --- otherwise concen. remains the same | ||
| 242 |
2/2✓ Branch 0 taken 869053 times.
✓ Branch 1 taken 348510 times.
|
1217563 | else for (p = 0; p < Nobjects[POLLUT]; p++) |
| 243 | { | ||
| 244 |
2/2✓ Branch 0 taken 88165 times.
✓ Branch 1 taken 780888 times.
|
869053 | if (Node[j].newDepth > ZeroDepth) |
| 245 | 88165 | Node[j].newQual[p] = Node[j].oldQual[p]; | |
| 246 | else | ||
| 247 | 780888 | Node[j].newQual[p] = 0.0; | |
| 248 | } | ||
| 249 | 2878459 | } | |
| 250 | |||
| 251 | //============================================================================= | ||
| 252 | |||
| 253 | 2758617 | void findLinkQual(int i, double tStep) | |
| 254 | // | ||
| 255 | // Input: i = link index | ||
| 256 | // tStep = routing time step (sec) | ||
| 257 | // Output: none | ||
| 258 | // Purpose: finds new quality in a link at end of the current time step. | ||
| 259 | // | ||
| 260 | { | ||
| 261 | int j, // upstream node index | ||
| 262 | k, // conduit index | ||
| 263 | p; // pollutant index | ||
| 264 | double wIn, // pollutant mass inflow rate (mass/sec) | ||
| 265 | qIn, // inflow rate (cfs) | ||
| 266 | qSeep, // rate of seepage loss (cfs) | ||
| 267 | v1, // link volume at start of time step (ft3) | ||
| 268 | v2, // link volume at end of time step (ft3) | ||
| 269 | c1, // current concentration within link (mass/ft3) | ||
| 270 | c2, // new concentration within link (mass/ft3) | ||
| 271 | vEvap, // volume lost to evaporation (ft3) | ||
| 272 | vLosses, // evap. + seepage volume loss (ft3) | ||
| 273 | fEvap, // evaporation concentration factor | ||
| 274 | barrels; // number of barrels in conduit | ||
| 275 | |||
| 276 | // --- identify index of upstream node | ||
| 277 | 2758617 | j = Link[i].node1; | |
| 278 |
2/2✓ Branch 0 taken 9817 times.
✓ Branch 1 taken 2748800 times.
|
2758617 | if ( Link[i].newFlow < 0.0 ) j = Link[i].node2; |
| 279 | |||
| 280 | // --- link quality is that of upstream node when | ||
| 281 | // link is not a conduit or is a dummy link | ||
| 282 |
3/4✓ Branch 0 taken 2679993 times.
✓ Branch 1 taken 78624 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2679993 times.
|
2758617 | if ( Link[i].type != CONDUIT || Link[i].xsect.type == DUMMY ) |
| 283 | { | ||
| 284 |
2/2✓ Branch 0 taken 206490 times.
✓ Branch 1 taken 78624 times.
|
285114 | for (p = 0; p < Nobjects[POLLUT]; p++) |
| 285 | { | ||
| 286 | 206490 | Link[i].newQual[p] = Node[j].newQual[p]; | |
| 287 | } | ||
| 288 | 78624 | return; | |
| 289 | } | ||
| 290 | |||
| 291 | // --- get flow rates and evaporation loss | ||
| 292 | 2679993 | k = Link[i].subIndex; | |
| 293 | 2679993 | barrels = Conduit[k].barrels; | |
| 294 | 2679993 | qIn = fabs(Conduit[k].q1) * barrels; | |
| 295 | 2679993 | qSeep = Conduit[k].seepLossRate * barrels; | |
| 296 | 2679993 | vEvap = Conduit[k].evapLossRate * barrels * tStep; | |
| 297 | |||
| 298 | // --- get starting and ending volumes | ||
| 299 | 2679993 | v1 = Link[i].oldVolume; | |
| 300 | 2679993 | v2 = Link[i].newVolume; | |
| 301 | 2679993 | vLosses = qSeep*tStep + vEvap; | |
| 302 | |||
| 303 | // --- compute factor by which concentrations are increased due to | ||
| 304 | // evaporation loss | ||
| 305 | 2679993 | fEvap = 1.0; | |
| 306 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2679993 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2679993 | if ( vEvap > 0.0 && v1 > ZeroVolume ) fEvap += vEvap / v1; |
| 307 | |||
| 308 | // --- Steady Flow routing requires special treatment | ||
| 309 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2679993 times.
|
2679993 | if ( RouteModel == SF ) |
| 310 | { | ||
| 311 | ✗ | findSFLinkQual(i, qSeep, fEvap, tStep); | |
| 312 | ✗ | return; | |
| 313 | } | ||
| 314 | |||
| 315 | // --- adjust inflow to compensate for volume change under Dynamic | ||
| 316 | // Wave routing (which produces just a single (out)flow rate | ||
| 317 | // for a conduit) | ||
| 318 |
1/2✓ Branch 0 taken 2679993 times.
✗ Branch 1 not taken.
|
2679993 | if ( RouteModel == DW ) |
| 319 | { | ||
| 320 | 2679993 | qIn = qIn + (v2 + vLosses - v1) / tStep; | |
| 321 |
2/2✓ Branch 0 taken 2489658 times.
✓ Branch 1 taken 190335 times.
|
2679993 | qIn = MAX(qIn, 0.0); |
| 322 | } | ||
| 323 | |||
| 324 | // --- examine each pollutant | ||
| 325 |
2/2✓ Branch 0 taken 6042941 times.
✓ Branch 1 taken 2679993 times.
|
8722934 | for (p = 0; p < Nobjects[POLLUT]; p++) |
| 326 | { | ||
| 327 | // --- start with concen. at start of time step | ||
| 328 | 6042941 | c1 = Link[i].oldQual[p]; | |
| 329 | |||
| 330 | // --- update mass balance accounting for seepage loss | ||
| 331 | 6042941 | massbal_addSeepageLoss(p, qSeep*c1); | |
| 332 | |||
| 333 | // --- increase concen. by evaporation factor | ||
| 334 | 6042941 | c1 *= fEvap; | |
| 335 | |||
| 336 | // --- reduce concen. by 1st-order reaction | ||
| 337 | 6042941 | c2 = getReactedQual(p, c1, v1, tStep); | |
| 338 | |||
| 339 | // --- mix resulting contents with inflow from upstream node | ||
| 340 | 6042941 | wIn = Node[j].newQual[p]*qIn; | |
| 341 | 6042941 | c2 = getMixedQual(c2, v1, wIn, qIn, tStep); | |
| 342 | |||
| 343 | // --- set concen. to zero if remaining volume is negligible | ||
| 344 |
4/4✓ Branch 0 taken 5103218 times.
✓ Branch 1 taken 939723 times.
✓ Branch 2 taken 915247 times.
✓ Branch 3 taken 4187971 times.
|
6042941 | if ( v2 < ZeroVolume || Link[i].newDepth <= ZeroDepth) |
| 345 | { | ||
| 346 | 1854970 | massbal_addToFinalStorage(p, c2 * v2); | |
| 347 | 1854970 | c2 = 0.0; | |
| 348 | } | ||
| 349 | |||
| 350 | // --- assign new concen. to link | ||
| 351 | 6042941 | Link[i].newQual[p] = c2; | |
| 352 | } | ||
| 353 | } | ||
| 354 | |||
| 355 | //============================================================================= | ||
| 356 | |||
| 357 | ✗ | void findSFLinkQual(int i, double qSeep, double fEvap, double tStep) | |
| 358 | // | ||
| 359 | // Input: i = link index | ||
| 360 | // tStep = routing time step (sec) | ||
| 361 | // Output: none | ||
| 362 | // Purpose: finds new quality in a link at end of the current time step for | ||
| 363 | // Steady Flow routing. | ||
| 364 | // | ||
| 365 | { | ||
| 366 | ✗ | int j = Link[i].node1; | |
| 367 | int p; | ||
| 368 | double c1, c2; | ||
| 369 | double lossRate; | ||
| 370 | |||
| 371 | // --- examine each pollutant | ||
| 372 | ✗ | for (p = 0; p < Nobjects[POLLUT]; p++) | |
| 373 | { | ||
| 374 | // --- conduit's quality equals upstream node quality | ||
| 375 | ✗ | c1 = Node[j].newQual[p]; | |
| 376 | |||
| 377 | // --- update mass balance accounting for seepage loss | ||
| 378 | ✗ | massbal_addSeepageLoss(p, qSeep*c1); | |
| 379 | |||
| 380 | // --- increase concen. by evaporation factor | ||
| 381 | ✗ | c1 *= fEvap; | |
| 382 | |||
| 383 | // --- apply first-order decay over travel time | ||
| 384 | ✗ | c2 = c1; | |
| 385 | ✗ | if ( Pollut[p].kDecay > 0.0 ) | |
| 386 | { | ||
| 387 | ✗ | c2 = c1 * exp(-Pollut[p].kDecay * tStep); | |
| 388 | ✗ | c2 = MAX(0.0, c2); | |
| 389 | ✗ | lossRate = (c1 - c2) * Link[i].newFlow; | |
| 390 | ✗ | massbal_addReactedMass(p, lossRate); | |
| 391 | } | ||
| 392 | ✗ | Link[i].newQual[p] = c2; | |
| 393 | } | ||
| 394 | ✗ | } | |
| 395 | |||
| 396 | //============================================================================= | ||
| 397 | |||
| 398 | 24147 | void findStorageQual(int j, double tStep) | |
| 399 | // | ||
| 400 | // Input: j = node index | ||
| 401 | // tStep = routing time step (sec) | ||
| 402 | // Output: none | ||
| 403 | // Purpose: finds new quality in a node with storage volume. | ||
| 404 | // | ||
| 405 | { | ||
| 406 | int p, // pollutant index | ||
| 407 | k; // storage unit index | ||
| 408 | double qIn, // inflow rate (cfs) | ||
| 409 | wIn, // pollutant mass inflow rate (mass) | ||
| 410 | v1, // volume at start of time step (ft3) | ||
| 411 | c1, // initial pollutant concentration (mass/ft3) | ||
| 412 | c2, // final pollutant concentration (mass/ft3) | ||
| 413 | 24147 | qExfil = 0.0, // exfiltration rate from storage unit (cfs) | |
| 414 | 24147 | vEvap = 0.0, // evaporation loss from storage unit (ft3) | |
| 415 | 24147 | fEvap = 1.0; // evaporation concentration factor | |
| 416 | |||
| 417 | // --- get inflow rate & initial volume | ||
| 418 | 24147 | qIn = Node[j].qualInflow; | |
| 419 | 24147 | v1 = Node[j].oldVolume; | |
| 420 | |||
| 421 | // -- for storage nodes | ||
| 422 |
1/2✓ Branch 0 taken 24147 times.
✗ Branch 1 not taken.
|
24147 | if ( Node[j].type == STORAGE ) |
| 423 | { | ||
| 424 | // --- update hydraulic residence time | ||
| 425 | // (HRT can be used in treatment functions) | ||
| 426 | 24147 | updateHRT(j, Node[j].oldVolume, qIn, tStep); | |
| 427 | |||
| 428 | // --- get exfiltration rate and evaporation loss | ||
| 429 | 24147 | k = Node[j].subIndex; | |
| 430 | 24147 | qExfil = Storage[k].exfilLoss / tStep; | |
| 431 | 24147 | vEvap = Storage[k].evapLoss; | |
| 432 | |||
| 433 | // --- compute factor by which concentrations are increased due to | ||
| 434 | // evaporation loss (avoiding huge factors as storage unit | ||
| 435 | // dries out completely) | ||
| 436 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 24147 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
24147 | if ( vEvap > 0.0 && v1 > ZeroVolume ) fEvap += vEvap / v1; |
| 437 | } | ||
| 438 | |||
| 439 | // --- for each pollutant | ||
| 440 |
2/2✓ Branch 0 taken 43667 times.
✓ Branch 1 taken 24147 times.
|
67814 | for (p = 0; p < Nobjects[POLLUT]; p++) |
| 441 | { | ||
| 442 | // --- start with concen. at start of time step | ||
| 443 | 43667 | c1 = Node[j].oldQual[p]; | |
| 444 | |||
| 445 | // --- update mass balance accounting for exfiltration loss | ||
| 446 | 43667 | massbal_addSeepageLoss(p, qExfil*c1); | |
| 447 | |||
| 448 | // --- increase concen. by evaporation factor | ||
| 449 | 43667 | c1 *= fEvap; | |
| 450 | |||
| 451 | // --- apply first order reaction only if no separate treatment function | ||
| 452 |
2/2✓ Branch 0 taken 5410 times.
✓ Branch 1 taken 38257 times.
|
43667 | if ( Node[j].treatment == NULL || |
| 453 |
2/2✓ Branch 0 taken 3787 times.
✓ Branch 1 taken 1623 times.
|
5410 | Node[j].treatment[p].equation == NULL ) |
| 454 | { | ||
| 455 | 42044 | c1 = getReactedQual(p, c1, v1, tStep); | |
| 456 | } | ||
| 457 | |||
| 458 | // --- mix resulting contents with inflow from all sources | ||
| 459 | // (temporarily accumulated in Node[j].newQual) | ||
| 460 | 43667 | wIn = Node[j].newQual[p]; | |
| 461 | 43667 | c2 = getMixedQual(c1, v1, wIn, qIn, tStep); | |
| 462 | |||
| 463 | // --- set concen. to zero if remaining volume & inflow is negligible | ||
| 464 |
2/2✓ Branch 0 taken 39233 times.
✓ Branch 1 taken 4434 times.
|
43667 | if ((Node[j].newVolume <= ZeroVolume || |
| 465 |
4/4✓ Branch 0 taken 388 times.
✓ Branch 1 taken 38845 times.
✓ Branch 2 taken 4006 times.
✓ Branch 3 taken 816 times.
|
43667 | Node[j].newDepth <= ZeroDepth) && qIn <= ZERO) |
| 466 | { | ||
| 467 | 4006 | massbal_addToFinalStorage(p, c2 * Node[j].newVolume); | |
| 468 | 4006 | c2 = 0.0; | |
| 469 | } | ||
| 470 | |||
| 471 | // --- assign new concen. to node | ||
| 472 | 43667 | Node[j].newQual[p] = c2; | |
| 473 | } | ||
| 474 | 24147 | } | |
| 475 | |||
| 476 | //============================================================================= | ||
| 477 | |||
| 478 | 24147 | void updateHRT(int j, double v, double q, double tStep) | |
| 479 | // | ||
| 480 | // Input: j = node index | ||
| 481 | // v = storage volume (ft3) | ||
| 482 | // q = inflow rate (cfs) | ||
| 483 | // tStep = time step (sec) | ||
| 484 | // Output: none | ||
| 485 | // Purpose: updates hydraulic residence time (i.e., water age) at a | ||
| 486 | // storage node. | ||
| 487 | // | ||
| 488 | { | ||
| 489 | 24147 | int k = Node[j].subIndex; | |
| 490 | 24147 | double hrt = Storage[k].hrt; | |
| 491 |
2/2✓ Branch 0 taken 1009 times.
✓ Branch 1 taken 23138 times.
|
24147 | if ( v < ZERO ) hrt = 0.0; |
| 492 | 23138 | else hrt = (hrt + tStep) * v / (v + q*tStep); | |
| 493 |
1/2✓ Branch 0 taken 24147 times.
✗ Branch 1 not taken.
|
24147 | Storage[k].hrt = MAX(hrt, 0.0); |
| 494 | 24147 | } | |
| 495 | |||
| 496 | //============================================================================= | ||
| 497 | |||
| 498 | 6084985 | double getReactedQual(int p, double c, double v1, double tStep) | |
| 499 | // | ||
| 500 | // Input: p = pollutant index | ||
| 501 | // c = initial concentration (mass/ft3) | ||
| 502 | // v1 = initial volume (ft3) | ||
| 503 | // tStep = time step (sec) | ||
| 504 | // Output: none | ||
| 505 | // Purpose: applies a first order reaction to a pollutant over a given | ||
| 506 | // time step. | ||
| 507 | // | ||
| 508 | { | ||
| 509 | double c2, lossRate; | ||
| 510 | 6084985 | double kDecay = Pollut[p].kDecay; | |
| 511 | |||
| 512 |
1/2✓ Branch 0 taken 6084985 times.
✗ Branch 1 not taken.
|
6084985 | if ( kDecay == 0.0 ) return c; |
| 513 | ✗ | c2 = c * (1.0 - kDecay * tStep); | |
| 514 | ✗ | c2 = MAX(0.0, c2); | |
| 515 | ✗ | lossRate = (c - c2) * v1 / tStep; | |
| 516 | ✗ | massbal_addReactedMass(p, lossRate); | |
| 517 | ✗ | return c2; | |
| 518 | } | ||
| 519 | |||
| 520 |