infil.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // infil.c | ||
| 3 | // | ||
| 4 | // Project: EPA SWMM5 | ||
| 5 | // Version: 5.2 | ||
| 6 | // Date: 11/01/21 (Build 5.2.0) | ||
| 7 | // Author: L. Rossman | ||
| 8 | // | ||
| 9 | // Infiltration functions. | ||
| 10 | // | ||
| 11 | // Update History | ||
| 12 | // ============== | ||
| 13 | // Build 5.1.007: | ||
| 14 | // - Revised formula for infiltration capacity recovery for the Modified | ||
| 15 | // Horton method. | ||
| 16 | // - The Green-Ampt functions were re-written. | ||
| 17 | // Build 5.1.008: | ||
| 18 | // - Monthly adjustment factors applied to hydraulic conductivity. | ||
| 19 | // Build 5.1.010: | ||
| 20 | // - Support for Modified Green Ampt model added. | ||
| 21 | // - Green-Ampt initial recovery time set to 0. | ||
| 22 | // Build 5.1.011: | ||
| 23 | // - Monthly hydraulic conductivity factor also applied to Fu parameter | ||
| 24 | // for Green-Ampt infiltration. | ||
| 25 | // - Prevented computed Horton infiltration from dropping below 0. | ||
| 26 | // Build 5.1.013: | ||
| 27 | // - Support added for subcatchment-specific time patterns that adjust | ||
| 28 | // hydraulic conductivity. | ||
| 29 | // Build 5.1.015: | ||
| 30 | // - Support added for multiple infiltration methods within a project. | ||
| 31 | // Build 5.2.0: | ||
| 32 | // - Additional validity check for G-A initial deficit added. | ||
| 33 | // - New error message 235 added for invalid infiltration parameters. | ||
| 34 | // - Conversion of runon to ponded depth fixed for Curve Number infiltration. | ||
| 35 | //----------------------------------------------------------------------------- | ||
| 36 | #define _CRT_SECURE_NO_DEPRECATE | ||
| 37 | |||
| 38 | #include <math.h> | ||
| 39 | #include <stdlib.h> | ||
| 40 | #include "headers.h" | ||
| 41 | #include "infil.h" | ||
| 42 | |||
| 43 | //----------------------------------------------------------------------------- | ||
| 44 | // Local Variables | ||
| 45 | //----------------------------------------------------------------------------- | ||
| 46 | typedef union TInfil { | ||
| 47 | THorton horton; | ||
| 48 | TGrnAmpt grnAmpt; | ||
| 49 | TCurveNum curveNum; | ||
| 50 | } TInfil; | ||
| 51 | TInfil *Infil; | ||
| 52 | |||
| 53 | static double Fumax; // saturated water volume in upper soil zone (ft) | ||
| 54 | static double InfilFactor; | ||
| 55 | |||
| 56 | //----------------------------------------------------------------------------- | ||
| 57 | // External Functions (declared in infil.h) | ||
| 58 | //----------------------------------------------------------------------------- | ||
| 59 | // infil_create (called by createObjects in project.c) | ||
| 60 | // infil_delete (called by deleteObjects in project.c) | ||
| 61 | // infil_readParams (called by input_readLine) | ||
| 62 | // infil_initState (called by subcatch_initState) | ||
| 63 | // infil_getState (called by writeRunoffFile in hotstart.c) | ||
| 64 | // infil_setState (called by readRunoffFile in hotstart.c) | ||
| 65 | // infil_getInfil (called by getSubareaRunoff in subcatch.c) | ||
| 66 | |||
| 67 | // Called locally and by storage node methods in node.c | ||
| 68 | // grnampt_setParams | ||
| 69 | // grnampt_initState | ||
| 70 | // grnampt_getInfil | ||
| 71 | |||
| 72 | //----------------------------------------------------------------------------- | ||
| 73 | // Local functions | ||
| 74 | //----------------------------------------------------------------------------- | ||
| 75 | static int horton_setParams(THorton *infil, double p[]); | ||
| 76 | static void horton_initState(THorton *infil); | ||
| 77 | static void horton_getState(THorton *infil, double x[]); | ||
| 78 | static void horton_setState(THorton *infil, double x[]); | ||
| 79 | static double horton_getInfil(THorton *infil, double tstep, double irate, | ||
| 80 | double depth); | ||
| 81 | static double modHorton_getInfil(THorton *infil, double tstep, double irate, | ||
| 82 | double depth); | ||
| 83 | |||
| 84 | static void grnampt_getState(TGrnAmpt *infil, double x[]); | ||
| 85 | static void grnampt_setState(TGrnAmpt *infil, double x[]); | ||
| 86 | static double grnampt_getUnsatInfil(TGrnAmpt *infil, double tstep, | ||
| 87 | double irate, double depth, int modelType); | ||
| 88 | static double grnampt_getSatInfil(TGrnAmpt *infil, double tstep, | ||
| 89 | double irate, double depth); | ||
| 90 | static double grnampt_getF2(double f1, double c1, double ks, double ts); | ||
| 91 | |||
| 92 | static int curvenum_setParams(TCurveNum *infil, double p[]); | ||
| 93 | static void curvenum_initState(TCurveNum *infil); | ||
| 94 | static void curvenum_getState(TCurveNum *infil, double x[]); | ||
| 95 | static void curvenum_setState(TCurveNum *infil, double x[]); | ||
| 96 | static double curvenum_getInfil(TCurveNum *infil, double tstep, double irate, | ||
| 97 | double depth); | ||
| 98 | |||
| 99 | //============================================================================= | ||
| 100 | |||
| 101 | 58 | void infil_create(int n) | |
| 102 | // | ||
| 103 | // Purpose: creates an array of infiltration objects. | ||
| 104 | // Input: n = number of subcatchments | ||
| 105 | // Output: none | ||
| 106 | // | ||
| 107 | { | ||
| 108 | 58 | Infil = (TInfil *) calloc(n, sizeof(TInfil)); | |
| 109 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 58 times.
|
58 | if (Infil == NULL) ErrorCode = ERR_MEMORY; |
| 110 | 58 | InfilFactor = 1.0; | |
| 111 | 58 | return; | |
| 112 | } | ||
| 113 | |||
| 114 | //============================================================================= | ||
| 115 | |||
| 116 | 58 | void infil_delete() | |
| 117 | // | ||
| 118 | // Purpose: deletes infiltration objects associated with subcatchments | ||
| 119 | // Input: none | ||
| 120 | // Output: none | ||
| 121 | // | ||
| 122 | { | ||
| 123 |
1/2✓ Branch 0 taken 58 times.
✗ Branch 1 not taken.
|
58 | FREE(Infil); |
| 124 | 58 | } | |
| 125 | |||
| 126 | //============================================================================= | ||
| 127 | |||
| 128 | 2389 | int infil_readParams(int m, char* tok[], int ntoks) | |
| 129 | // | ||
| 130 | // Input: m = default infiltration model | ||
| 131 | // tok[] = array of string tokens | ||
| 132 | // ntoks = number of tokens | ||
| 133 | // Output: returns an error code | ||
| 134 | // Purpose: sets infiltration parameters from a line of input data. | ||
| 135 | // | ||
| 136 | // Format of data line is: | ||
| 137 | // subcatch p1 p2 ... (infilMethod) | ||
| 138 | { | ||
| 139 | int i, j, n, status; | ||
| 140 | double x[5]; | ||
| 141 | |||
| 142 | // --- check that subcatchment exists | ||
| 143 | 2389 | j = project_findObject(SUBCATCH, tok[0]); | |
| 144 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2389 times.
|
2389 | if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]); |
| 145 | |||
| 146 | // --- check for infiltration method keyword is last token | ||
| 147 | 2389 | i = findmatch(tok[ntoks-1], InfilModelWords); | |
| 148 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2389 times.
|
2389 | if ( i >= 0 ) |
| 149 | { | ||
| 150 | ✗ | m = i; | |
| 151 | ✗ | --ntoks; | |
| 152 | } | ||
| 153 | |||
| 154 | // --- number of input tokens depends on infiltration model m | ||
| 155 |
2/2✓ Branch 0 taken 59 times.
✓ Branch 1 taken 2330 times.
|
2389 | if ( m == HORTON ) n = 5; |
| 156 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2329 times.
|
2330 | else if ( m == MOD_HORTON ) n = 5; |
| 157 |
2/2✓ Branch 0 taken 2319 times.
✓ Branch 1 taken 10 times.
|
2329 | else if ( m == GREEN_AMPT ) n = 4; |
| 158 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9 times.
|
10 | else if ( m == MOD_GREEN_AMPT ) n = 4; |
| 159 |
1/2✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
|
9 | else if ( m == CURVE_NUMBER ) n = 4; |
| 160 | ✗ | else return 0; | |
| 161 | |||
| 162 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2389 times.
|
2389 | if ( ntoks < n ) return error_setInpError(ERR_ITEMS, ""); |
| 163 | |||
| 164 | // --- parse numerical values from tokens | ||
| 165 |
2/2✓ Branch 0 taken 11945 times.
✓ Branch 1 taken 2389 times.
|
14334 | for (i = 0; i < 5; i++) x[i] = 0.0; |
| 166 |
2/2✓ Branch 0 taken 7227 times.
✓ Branch 1 taken 2389 times.
|
9616 | for (i = 1; i < n; i++) |
| 167 | { | ||
| 168 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 7227 times.
|
7227 | if (!getDouble(tok[i], &x[i - 1])) |
| 169 | ✗ | return error_setInpError(ERR_NUMBER, tok[i]); | |
| 170 | } | ||
| 171 | |||
| 172 | // --- special case for Horton infil. - last parameter is optional | ||
| 173 |
5/6✓ Branch 0 taken 2330 times.
✓ Branch 1 taken 59 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 2329 times.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
|
2389 | if ( (m == HORTON || m == MOD_HORTON) && ntoks > n ) |
| 174 | { | ||
| 175 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
|
60 | if ( ! getDouble(tok[n], &x[n-1]) ) |
| 176 | ✗ | return error_setInpError(ERR_NUMBER, tok[n]); | |
| 177 | } | ||
| 178 | |||
| 179 | // --- assign parameter values to infil, infilModel object | ||
| 180 | 2389 | Subcatch[j].infil = j; | |
| 181 | 2389 | Subcatch[j].infilModel = m; | |
| 182 |
3/4✓ Branch 0 taken 60 times.
✓ Branch 1 taken 2320 times.
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
|
2389 | switch (m) |
| 183 | { | ||
| 184 | 60 | case HORTON: | |
| 185 | 60 | case MOD_HORTON: status = horton_setParams(&Infil[j].horton, x); | |
| 186 | 60 | break; | |
| 187 | 2320 | case GREEN_AMPT: | |
| 188 | case MOD_GREEN_AMPT: | ||
| 189 | 2320 | status = grnampt_setParams(&Infil[j].grnAmpt, x); | |
| 190 | 2320 | break; | |
| 191 | 9 | case CURVE_NUMBER: status = curvenum_setParams(&Infil[j].curveNum, x); | |
| 192 | 9 | break; | |
| 193 | ✗ | default: status = TRUE; | |
| 194 | } | ||
| 195 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2389 times.
|
2389 | if ( !status ) return error_setInpError(ERR_INFIL_PARAMS, ""); |
| 196 | 2389 | return 0; | |
| 197 | } | ||
| 198 | |||
| 199 | //============================================================================= | ||
| 200 | |||
| 201 | 2389 | void infil_initState(int j) | |
| 202 | // | ||
| 203 | // Input: j = subcatchment index | ||
| 204 | // Output: none | ||
| 205 | // Purpose: initializes state of infiltration for a subcatchment. | ||
| 206 | // | ||
| 207 | { | ||
| 208 |
3/4✓ Branch 0 taken 60 times.
✓ Branch 1 taken 2320 times.
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
|
2389 | switch (Subcatch[j].infilModel) |
| 209 | { | ||
| 210 | 60 | case HORTON: | |
| 211 | 60 | case MOD_HORTON: horton_initState(&Infil[j].horton); break; | |
| 212 | 2320 | case GREEN_AMPT: | |
| 213 | case MOD_GREEN_AMPT: | ||
| 214 | 2320 | grnampt_initState(&Infil[j].grnAmpt); break; | |
| 215 | 9 | case CURVE_NUMBER: curvenum_initState(&Infil[j].curveNum); break; | |
| 216 | } | ||
| 217 | 2389 | } | |
| 218 | |||
| 219 | //============================================================================= | ||
| 220 | |||
| 221 | 329 | void infil_getState(int j, double x[]) | |
| 222 | // | ||
| 223 | // Input: j = subcatchment index | ||
| 224 | // Output: x = subcatchment's infiltration state | ||
| 225 | // Purpose: retrieves the current infiltration state for a subcatchment. | ||
| 226 | // | ||
| 227 | { | ||
| 228 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 329 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
329 | switch (Subcatch[j].infilModel) |
| 229 | { | ||
| 230 | ✗ | case HORTON: | |
| 231 | ✗ | case MOD_HORTON: horton_getState(&Infil[j].horton, x); break; | |
| 232 | 329 | case GREEN_AMPT: | |
| 233 | case MOD_GREEN_AMPT: | ||
| 234 | 329 | grnampt_getState(&Infil[j].grnAmpt, x); break; | |
| 235 | ✗ | case CURVE_NUMBER: curvenum_getState(&Infil[j].curveNum, x); break; | |
| 236 | } | ||
| 237 | 329 | } | |
| 238 | |||
| 239 | //============================================================================= | ||
| 240 | |||
| 241 | 329 | void infil_setState(int j, double x[]) | |
| 242 | // | ||
| 243 | // Input: j = subcatchment index | ||
| 244 | // m = infiltration method code | ||
| 245 | // Output: none | ||
| 246 | // Purpose: sets the current infiltration state for a subcatchment. | ||
| 247 | // | ||
| 248 | { | ||
| 249 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 329 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
329 | switch (Subcatch[j].infilModel) |
| 250 | { | ||
| 251 | ✗ | case HORTON: | |
| 252 | ✗ | case MOD_HORTON: horton_setState(&Infil[j].horton, x); break; | |
| 253 | 329 | case GREEN_AMPT: | |
| 254 | case MOD_GREEN_AMPT: | ||
| 255 | 329 | grnampt_setState(&Infil[j].grnAmpt, x); break; | |
| 256 | ✗ | case CURVE_NUMBER: curvenum_setState(&Infil[j].curveNum, x); break; | |
| 257 | } | ||
| 258 | 329 | } | |
| 259 | |||
| 260 | //============================================================================= | ||
| 261 | |||
| 262 | 4296115 | void infil_setInfilFactor(int j) | |
| 263 | // | ||
| 264 | // Input: j = subcatchment index | ||
| 265 | // Output: none | ||
| 266 | // Purpose: assigns a value to the infiltration adjustment factor. | ||
| 267 | { | ||
| 268 | int m; | ||
| 269 | int p; | ||
| 270 | |||
| 271 | // ... set factor to the global conductivity adjustment factor | ||
| 272 | 4296115 | InfilFactor = Adjust.hydconFactor; | |
| 273 | |||
| 274 | // ... override global factor with subcatchment's adjustment if assigned | ||
| 275 |
2/2✓ Branch 0 taken 3229694 times.
✓ Branch 1 taken 1066421 times.
|
4296115 | if (j >= 0) |
| 276 | { | ||
| 277 | 3229694 | p = Subcatch[j].infilPattern; | |
| 278 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 3229694 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
3229694 | if (p >= 0 && Pattern[p].type == MONTHLY_PATTERN) |
| 279 | { | ||
| 280 | ✗ | m = datetime_monthOfYear(getDateTime(OldRunoffTime)) - 1; | |
| 281 | ✗ | InfilFactor = Pattern[p].factor[m]; | |
| 282 | } | ||
| 283 | } | ||
| 284 | 4296115 | } | |
| 285 | |||
| 286 | //============================================================================= | ||
| 287 | |||
| 288 | 3228131 | double infil_getInfil(int j, double tstep, double rainfall, | |
| 289 | double runon, double depth) | ||
| 290 | // | ||
| 291 | // Input: j = subcatchment index | ||
| 292 | // tstep = runoff time step (sec) | ||
| 293 | // rainfall = rainfall rate (ft/sec) | ||
| 294 | // runon = runon rate from other sub-areas or subcatchments (ft/sec) | ||
| 295 | // depth = depth of surface water on subcatchment (ft) | ||
| 296 | // Output: returns infiltration rate (ft/sec) | ||
| 297 | // Purpose: computes infiltration rate depending on infiltration method. | ||
| 298 | // | ||
| 299 | { | ||
| 300 |
4/5✓ Branch 0 taken 88921 times.
✓ Branch 1 taken 121 times.
✓ Branch 2 taken 3138162 times.
✓ Branch 3 taken 927 times.
✗ Branch 4 not taken.
|
3228131 | switch (Subcatch[j].infilModel) |
| 301 | { | ||
| 302 | 88921 | case HORTON: | |
| 303 | 88921 | return horton_getInfil(&Infil[j].horton, tstep, rainfall+runon, depth); | |
| 304 | |||
| 305 | 121 | case MOD_HORTON: | |
| 306 | 121 | return modHorton_getInfil(&Infil[j].horton, tstep, rainfall+runon, | |
| 307 | depth); | ||
| 308 | |||
| 309 | 3138162 | case GREEN_AMPT: | |
| 310 | case MOD_GREEN_AMPT: | ||
| 311 | 3138162 | return grnampt_getInfil(&Infil[j].grnAmpt, tstep, rainfall+runon, depth, | |
| 312 | 3138162 | Subcatch[j].infilModel); | |
| 313 | |||
| 314 | 927 | case CURVE_NUMBER: | |
| 315 | 927 | depth += runon * tstep; | |
| 316 | 927 | return curvenum_getInfil(&Infil[j].curveNum, tstep, rainfall, depth); | |
| 317 | |||
| 318 | ✗ | default: | |
| 319 | ✗ | return 0.0; | |
| 320 | } | ||
| 321 | } | ||
| 322 | |||
| 323 | //============================================================================= | ||
| 324 | |||
| 325 | 60 | int horton_setParams(THorton *infil, double p[]) | |
| 326 | // | ||
| 327 | // Input: infil = ptr. to Horton infiltration object | ||
| 328 | // p[] = array of parameter values | ||
| 329 | // Output: returns TRUE if parameters are valid, FALSE otherwise | ||
| 330 | // Purpose: assigns Horton infiltration parameters to a subcatchment. | ||
| 331 | // | ||
| 332 | { | ||
| 333 | int k; | ||
| 334 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 300 times.
✓ Branch 2 taken 300 times.
✓ Branch 3 taken 60 times.
|
360 | for (k = 0; k < 5; k++) if ( p[k] < 0.0 ) return FALSE; |
| 335 | |||
| 336 | // --- max. & min. infil rates (ft/sec) | ||
| 337 | 60 | infil->f0 = p[0] / UCF(RAINFALL); | |
| 338 | 60 | infil->fmin = p[1] / UCF(RAINFALL); | |
| 339 | |||
| 340 | // --- convert decay const. to 1/sec | ||
| 341 | 60 | infil->decay = p[2] / 3600.; | |
| 342 | |||
| 343 | // --- convert drying time (days) to a regeneration const. (1/sec) | ||
| 344 | // assuming that former is time to reach 98% dry along an | ||
| 345 | // exponential drying curve | ||
| 346 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
|
60 | if (p[3] == 0.0 ) p[3] = TINY; |
| 347 | 60 | infil->regen = -log(1.0-0.98) / p[3] / SECperDAY; | |
| 348 | |||
| 349 | // --- optional max. infil. capacity (ft) (p[4] = 0 if no value supplied) | ||
| 350 | 60 | infil->Fmax = p[4] / UCF(RAINDEPTH); | |
| 351 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
|
60 | if ( infil->f0 < infil->fmin ) return FALSE; |
| 352 | 60 | return TRUE; | |
| 353 | } | ||
| 354 | |||
| 355 | //============================================================================= | ||
| 356 | |||
| 357 | 60 | void horton_initState(THorton *infil) | |
| 358 | // | ||
| 359 | // Input: infil = ptr. to Horton infiltration object | ||
| 360 | // Output: none | ||
| 361 | // Purpose: initializes time on Horton infiltration curve for a subcatchment. | ||
| 362 | // | ||
| 363 | { | ||
| 364 | 60 | infil->tp = 0.0; | |
| 365 | 60 | infil->Fe = 0.0; | |
| 366 | 60 | } | |
| 367 | |||
| 368 | //============================================================================= | ||
| 369 | |||
| 370 | ✗ | void horton_getState(THorton *infil, double x[]) | |
| 371 | { | ||
| 372 | ✗ | x[0] = infil->tp; | |
| 373 | ✗ | x[1] = infil->Fe; | |
| 374 | ✗ | } | |
| 375 | |||
| 376 | ✗ | void horton_setState(THorton *infil, double x[]) | |
| 377 | { | ||
| 378 | ✗ | infil->tp = x[0]; | |
| 379 | ✗ | infil->Fe = x[1]; | |
| 380 | ✗ | } | |
| 381 | |||
| 382 | //============================================================================= | ||
| 383 | |||
| 384 | 88921 | double horton_getInfil(THorton *infil, double tstep, double irate, double depth) | |
| 385 | // | ||
| 386 | // Input: infil = ptr. to Horton infiltration object | ||
| 387 | // tstep = runoff time step (sec), | ||
| 388 | // irate = net "rainfall" rate (ft/sec), | ||
| 389 | // = rainfall + snowmelt + runon - evaporation | ||
| 390 | // depth = depth of ponded water (ft). | ||
| 391 | // Output: returns infiltration rate (ft/sec) | ||
| 392 | // Purpose: computes Horton infiltration for a subcatchment. | ||
| 393 | // | ||
| 394 | { | ||
| 395 | // --- assign local variables | ||
| 396 | int iter; | ||
| 397 | 88921 | double fa, fp = 0.0; | |
| 398 | double Fp, F1, t1, tlim, ex, kt; | ||
| 399 | double FF, FF1, r; | ||
| 400 | 88921 | double f0 = infil->f0 * InfilFactor; | |
| 401 | 88921 | double fmin = infil->fmin * InfilFactor; | |
| 402 | 88921 | double Fmax = infil->Fmax; | |
| 403 | 88921 | double tp = infil->tp; | |
| 404 | 88921 | double df = f0 - fmin; | |
| 405 | 88921 | double kd = infil->decay; | |
| 406 | 88921 | double kr = infil->regen * Evap.recoveryFactor; | |
| 407 | |||
| 408 | // --- special cases of no infil. or constant infil | ||
| 409 |
3/6✓ Branch 0 taken 88921 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 88921 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 88921 times.
|
88921 | if ( df < 0.0 || kd < 0.0 || kr < 0.0 ) return 0.0; |
| 410 |
3/4✓ Branch 0 taken 88825 times.
✓ Branch 1 taken 96 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 88825 times.
|
88921 | if ( df == 0.0 || kd == 0.0 ) |
| 411 | { | ||
| 412 | 96 | fp = f0; | |
| 413 | 96 | fa = irate + depth / tstep; | |
| 414 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 96 times.
|
96 | if ( fp > fa ) fp = fa; |
| 415 |
1/2✓ Branch 0 taken 96 times.
✗ Branch 1 not taken.
|
96 | return MAX(0.0, fp); |
| 416 | } | ||
| 417 | |||
| 418 | // --- compute water available for infiltration | ||
| 419 | 88825 | fa = irate + depth / tstep; | |
| 420 | |||
| 421 | // --- case where there is water to infiltrate | ||
| 422 |
2/2✓ Branch 0 taken 15550 times.
✓ Branch 1 taken 73275 times.
|
88825 | if ( fa > ZERO ) |
| 423 | { | ||
| 424 | // --- compute average infil. rate over time step | ||
| 425 | 15550 | t1 = tp + tstep; // future cumul. time | |
| 426 | 15550 | tlim = 16.0 / kd; // for tp >= tlim, f = fmin | |
| 427 |
2/2✓ Branch 0 taken 2435 times.
✓ Branch 1 taken 13115 times.
|
15550 | if ( tp >= tlim ) |
| 428 | { | ||
| 429 | 2435 | Fp = fmin * tp + df / kd; | |
| 430 | 2435 | F1 = Fp + fmin * tstep; | |
| 431 | } | ||
| 432 | else | ||
| 433 | { | ||
| 434 | 13115 | Fp = fmin * tp + df / kd * (1.0 - exp(-kd * tp)); | |
| 435 | 13115 | F1 = fmin * t1 + df / kd * (1.0 - exp(-kd * t1)); | |
| 436 | } | ||
| 437 | 15550 | fp = (F1 - Fp) / tstep; | |
| 438 |
2/2✓ Branch 0 taken 13842 times.
✓ Branch 1 taken 1708 times.
|
15550 | fp = MAX(fp, fmin); |
| 439 | |||
| 440 | // --- limit infil rate to available infil | ||
| 441 |
2/2✓ Branch 0 taken 13405 times.
✓ Branch 1 taken 2145 times.
|
15550 | if ( fp > fa ) fp = fa; |
| 442 | |||
| 443 | // --- if fp on flat portion of curve then increase tp by tstep | ||
| 444 |
2/2✓ Branch 0 taken 2477 times.
✓ Branch 1 taken 13073 times.
|
15550 | if ( t1 > tlim ) tp = t1; |
| 445 | |||
| 446 | // --- if infil < available capacity then increase tp by tstep | ||
| 447 |
2/2✓ Branch 0 taken 1541 times.
✓ Branch 1 taken 11532 times.
|
13073 | else if ( fp < fa ) tp = t1; |
| 448 | |||
| 449 | // --- if infil limited by available capacity then | ||
| 450 | // solve F(tp) - F1 = 0 using Newton-Raphson method | ||
| 451 | else | ||
| 452 | { | ||
| 453 | 11532 | F1 = Fp + fp * tstep; | |
| 454 | 11532 | tp = tp + tstep / 2.0; | |
| 455 |
1/2✓ Branch 0 taken 30161 times.
✗ Branch 1 not taken.
|
30161 | for ( iter=1; iter<=20; iter++ ) |
| 456 | { | ||
| 457 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 30161 times.
|
30161 | kt = MIN( 60.0, kd*tp ); |
| 458 | 30161 | ex = exp(-kt); | |
| 459 | 30161 | FF = fmin * tp + df / kd * (1.0 - ex) - F1; | |
| 460 | 30161 | FF1 = fmin + df * ex; | |
| 461 | 30161 | r = FF / FF1; | |
| 462 | 30161 | tp = tp - r; | |
| 463 |
2/2✓ Branch 0 taken 11532 times.
✓ Branch 1 taken 18629 times.
|
30161 | if ( fabs(r) <= 0.001 * tstep ) break; |
| 464 | } | ||
| 465 | } | ||
| 466 | |||
| 467 | // --- limit cumulative infiltration to Fmax | ||
| 468 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15550 times.
|
15550 | if ( Fmax > 0.0 ) |
| 469 | { | ||
| 470 | ✗ | if ( infil->Fe + fp * tstep > Fmax ) | |
| 471 | ✗ | fp = (Fmax - infil->Fe) / tstep; | |
| 472 | ✗ | fp = MAX(fp, 0.0); | |
| 473 | ✗ | infil->Fe += fp * tstep; | |
| 474 | } | ||
| 475 | } | ||
| 476 | |||
| 477 | // --- case where infil. capacity is regenerating; update tp. | ||
| 478 |
1/2✓ Branch 0 taken 73275 times.
✗ Branch 1 not taken.
|
73275 | else if (kr > 0.0) |
| 479 | { | ||
| 480 | 73275 | r = exp(-kr * tstep); | |
| 481 | 73275 | tp = 1.0 - exp(-kd * tp); | |
| 482 | 73275 | tp = -log(1.0 - r*tp) / kd; | |
| 483 | |||
| 484 | // reduction in cumulative infiltration | ||
| 485 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 73275 times.
|
73275 | if ( Fmax > 0.0 ) |
| 486 | { | ||
| 487 | ✗ | infil->Fe = fmin*tp + (df/kd)*(1.0 - exp(-kd*tp)); | |
| 488 | } | ||
| 489 | } | ||
| 490 | 88825 | infil->tp = tp; | |
| 491 | 88825 | return fp; | |
| 492 | } | ||
| 493 | |||
| 494 | //============================================================================= | ||
| 495 | |||
| 496 | 121 | double modHorton_getInfil(THorton *infil, double tstep, double irate, | |
| 497 | double depth) | ||
| 498 | // | ||
| 499 | // Input: infil = ptr. to Horton infiltration object | ||
| 500 | // tstep = runoff time step (sec), | ||
| 501 | // irate = net "rainfall" rate (ft/sec), | ||
| 502 | // = rainfall + snowmelt + runon | ||
| 503 | // depth = depth of ponded water (ft). | ||
| 504 | // Output: returns infiltration rate (ft/sec) | ||
| 505 | // Purpose: computes modified Horton infiltration for a subcatchment. | ||
| 506 | // | ||
| 507 | { | ||
| 508 | // --- assign local variables | ||
| 509 | 121 | double f = 0.0; | |
| 510 | double fp, fa; | ||
| 511 | 121 | double f0 = infil->f0 * InfilFactor; | |
| 512 | 121 | double fmin = infil->fmin * InfilFactor; | |
| 513 | 121 | double df = f0 - fmin; | |
| 514 | 121 | double kd = infil->decay; | |
| 515 | 121 | double kr = infil->regen * Evap.recoveryFactor; | |
| 516 | |||
| 517 | // --- special cases of no or constant infiltration | ||
| 518 |
3/6✓ Branch 0 taken 121 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 121 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 121 times.
|
121 | if ( df < 0.0 || kd < 0.0 || kr < 0.0 ) return 0.0; |
| 519 |
2/4✓ Branch 0 taken 121 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 121 times.
|
121 | if ( df == 0.0 || kd == 0.0 ) |
| 520 | { | ||
| 521 | ✗ | fp = f0; | |
| 522 | ✗ | fa = irate + depth / tstep; | |
| 523 | ✗ | if ( fp > fa ) fp = fa; | |
| 524 | ✗ | return MAX(0.0, fp); | |
| 525 | } | ||
| 526 | |||
| 527 | // --- compute water available for infiltration | ||
| 528 | 121 | fa = irate + depth / tstep; | |
| 529 | |||
| 530 | // --- case where there is water to infiltrate | ||
| 531 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 11 times.
|
121 | if ( fa > ZERO ) |
| 532 | { | ||
| 533 | // --- saturated condition | ||
| 534 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
110 | if ( infil->Fmax > 0.0 && infil->Fe >= infil->Fmax ) return 0.0; |
| 535 | |||
| 536 | // --- potential infiltration | ||
| 537 | 110 | fp = f0 - kd * infil->Fe; | |
| 538 |
1/2✓ Branch 0 taken 110 times.
✗ Branch 1 not taken.
|
110 | fp = MAX(fp, fmin); |
| 539 | |||
| 540 | // --- actual infiltration | ||
| 541 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 109 times.
|
110 | f = MIN(fa, fp); |
| 542 | |||
| 543 | // --- new cumulative infiltration minus seepage | ||
| 544 |
2/2✓ Branch 0 taken 109 times.
✓ Branch 1 taken 1 time.
|
110 | infil->Fe += MAX((f - fmin), 0.0) * tstep; |
| 545 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
110 | if ( infil->Fmax > 0.0 ) infil->Fe = MAX(infil->Fe, infil->Fmax); |
| 546 | } | ||
| 547 | |||
| 548 | // --- reduce cumulative infiltration for dry condition | ||
| 549 |
1/2✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
|
11 | else if (kr > 0.0) |
| 550 | { | ||
| 551 | 11 | infil->Fe *= exp(-kr * tstep); | |
| 552 |
1/2✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
|
11 | infil->Fe = MAX(infil->Fe, 0.0); |
| 553 | } | ||
| 554 | 121 | return f; | |
| 555 | } | ||
| 556 | |||
| 557 | //============================================================================= | ||
| 558 | |||
| 559 | ✗ | void grnampt_getParams(int j, double p[]) | |
| 560 | // | ||
| 561 | // Input: j = subcatchment index | ||
| 562 | // p[] = array of parameter values | ||
| 563 | // Output: none | ||
| 564 | // Purpose: retrieves Green-Ampt infiltration parameters for a subcatchment. | ||
| 565 | // | ||
| 566 | { | ||
| 567 | ✗ | p[0] = Infil[j].grnAmpt.S * UCF(RAINDEPTH); // Capillary suction head (ft) | |
| 568 | ✗ | p[1] = Infil[j].grnAmpt.Ks * UCF(RAINFALL); // Sat. hyd. conductivity (ft/sec) | |
| 569 | ✗ | p[2] = Infil[j].grnAmpt.IMDmax; // Max. init. moisture deficit | |
| 570 | ✗ | } | |
| 571 | |||
| 572 | //============================================================================= | ||
| 573 | |||
| 574 | 2334 | int grnampt_setParams(TGrnAmpt *infil, double p[]) | |
| 575 | // | ||
| 576 | // Input: infil = ptr. to Green-Ampt infiltration object | ||
| 577 | // p[] = array of parameter values | ||
| 578 | // Output: returns TRUE if parameters are valid, FALSE otherwise | ||
| 579 | // Purpose: assigns Green-Ampt infiltration parameters to a subcatchment. | ||
| 580 | // | ||
| 581 | { | ||
| 582 | double ksat; // sat. hyd. conductivity in in/hr | ||
| 583 | |||
| 584 |
4/8✓ Branch 0 taken 2334 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2334 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2334 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2334 times.
|
2334 | if ( p[0] < 0.0 || p[1] <= 0.0 || p[2] < 0.0 || p[2] > 1.0) return FALSE; |
| 585 | 2334 | infil->S = p[0] / UCF(RAINDEPTH); // Capillary suction head (ft) | |
| 586 | 2334 | infil->Ks = p[1] / UCF(RAINFALL); // Sat. hyd. conductivity (ft/sec) | |
| 587 | 2334 | infil->IMDmax = p[2]; // Max. init. moisture deficit | |
| 588 | |||
| 589 | // --- find depth of upper soil zone (ft) using Mein's eqn. | ||
| 590 | 2334 | ksat = infil->Ks * 12. * 3600.; | |
| 591 | 2334 | infil->Lu = 4.0 * sqrt(ksat) / 12.; | |
| 592 | 2334 | return TRUE; | |
| 593 | } | ||
| 594 | |||
| 595 | //============================================================================= | ||
| 596 | |||
| 597 | 2334 | void grnampt_initState(TGrnAmpt *infil) | |
| 598 | // | ||
| 599 | // Input: infil = ptr. to Green-Ampt infiltration object | ||
| 600 | // Output: none | ||
| 601 | // Purpose: initializes state of Green-Ampt infiltration for a subcatchment. | ||
| 602 | // | ||
| 603 | { | ||
| 604 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2334 times.
|
2334 | if (infil == NULL) return; |
| 605 | 2334 | infil->IMD = infil->IMDmax; | |
| 606 | 2334 | infil->Fu = 0.0; | |
| 607 | 2334 | infil->F = 0.0; | |
| 608 | 2334 | infil->Sat = FALSE; | |
| 609 | 2334 | infil->T = 0.0; | |
| 610 | } | ||
| 611 | |||
| 612 | 329 | void grnampt_getState(TGrnAmpt *infil, double x[]) | |
| 613 | { | ||
| 614 | 329 | x[0] = infil->IMD; | |
| 615 | 329 | x[1] = infil->F; | |
| 616 | 329 | x[2] = infil->Fu; | |
| 617 | 329 | x[3] = infil->Sat; | |
| 618 | 329 | x[4] = infil->T; | |
| 619 | 329 | } | |
| 620 | |||
| 621 | 329 | void grnampt_setState(TGrnAmpt *infil, double x[]) | |
| 622 | { | ||
| 623 | 329 | infil->IMD = x[0]; | |
| 624 | 329 | infil->F = x[1]; | |
| 625 | 329 | infil->Fu = x[2]; | |
| 626 | 329 | infil->Sat = (char)x[3]; | |
| 627 | 329 | infil->T = x[4]; | |
| 628 | 329 | } | |
| 629 | |||
| 630 | //============================================================================= | ||
| 631 | |||
| 632 | 3177235 | double grnampt_getInfil(TGrnAmpt *infil, double tstep, double irate, | |
| 633 | double depth, int modelType) | ||
| 634 | // | ||
| 635 | // Input: infil = ptr. to Green-Ampt infiltration object | ||
| 636 | // tstep = time step (sec), | ||
| 637 | // irate = net "rainfall" rate to upper zone (ft/sec); | ||
| 638 | // = rainfall + snowmelt + runon, | ||
| 639 | // does not include ponded water (added on below) | ||
| 640 | // depth = depth of ponded water (ft) | ||
| 641 | // modelType = either GREEN_AMPT or MOD_GREEN_AMPT | ||
| 642 | // Output: returns infiltration rate (ft/sec) | ||
| 643 | // Purpose: computes Green-Ampt infiltration for a subcatchment | ||
| 644 | // or a storage node. | ||
| 645 | // | ||
| 646 | { | ||
| 647 | // --- find saturated upper soil zone water volume | ||
| 648 | 3177235 | Fumax = infil->IMDmax * infil->Lu * sqrt(InfilFactor); | |
| 649 | |||
| 650 | // --- reduce time until next event | ||
| 651 | 3177235 | infil->T -= tstep; | |
| 652 | |||
| 653 | // --- use different procedures depending on upper soil zone saturation | ||
| 654 |
2/2✓ Branch 0 taken 61122 times.
✓ Branch 1 taken 3116113 times.
|
3177235 | if ( infil->Sat ) return grnampt_getSatInfil(infil, tstep, irate, depth); |
| 655 | 3116113 | else return grnampt_getUnsatInfil(infil, tstep, irate, depth, modelType); | |
| 656 | } | ||
| 657 | |||
| 658 | //============================================================================= | ||
| 659 | |||
| 660 | 3116113 | double grnampt_getUnsatInfil(TGrnAmpt *infil, double tstep, double irate, | |
| 661 | double depth, int modelType) | ||
| 662 | // | ||
| 663 | // Input: infil = ptr. to Green-Ampt infiltration object | ||
| 664 | // tstep = runoff time step (sec), | ||
| 665 | // irate = net "rainfall" rate to upper zone (ft/sec); | ||
| 666 | // = rainfall + snowmelt + runon, | ||
| 667 | // does not include ponded water (added on below) | ||
| 668 | // depth = depth of ponded water (ft) | ||
| 669 | // modelType = either GREEN_AMPT or MOD_GREEN_AMPT | ||
| 670 | // Output: returns infiltration rate (ft/sec) | ||
| 671 | // Purpose: computes Green-Ampt infiltration when upper soil zone is | ||
| 672 | // unsaturated. | ||
| 673 | // | ||
| 674 | { | ||
| 675 | double ia, c1, F2, dF, Fs, kr, ts; | ||
| 676 | 3116113 | double ks = infil->Ks * InfilFactor; | |
| 677 | 3116113 | double lu = infil->Lu * sqrt(InfilFactor); | |
| 678 | |||
| 679 | // --- get available infiltration rate (rainfall + ponded water) | ||
| 680 | 3116113 | ia = irate + depth / tstep; | |
| 681 |
2/2✓ Branch 0 taken 2302925 times.
✓ Branch 1 taken 813188 times.
|
3116113 | if ( ia < ZERO ) ia = 0.0; |
| 682 | |||
| 683 | // --- no rainfall so recover upper zone moisture | ||
| 684 |
2/2✓ Branch 0 taken 2302925 times.
✓ Branch 1 taken 813188 times.
|
3116113 | if ( ia == 0.0 ) |
| 685 | { | ||
| 686 |
2/2✓ Branch 0 taken 779365 times.
✓ Branch 1 taken 1523560 times.
|
2302925 | if ( infil->Fu <= 0.0 ) return 0.0; |
| 687 | 1523560 | kr = lu / 90000.0 * Evap.recoveryFactor; | |
| 688 | 1523560 | dF = kr * Fumax * tstep; | |
| 689 | 1523560 | infil->F -= dF; | |
| 690 | 1523560 | infil->Fu -= dF; | |
| 691 |
2/2✓ Branch 0 taken 2709 times.
✓ Branch 1 taken 1520851 times.
|
1523560 | if ( infil->Fu <= 0.0 ) |
| 692 | { | ||
| 693 | 2709 | infil->Fu = 0.0; | |
| 694 | 2709 | infil->F = 0.0; | |
| 695 | 2709 | infil->IMD = infil->IMDmax; | |
| 696 | 2709 | return 0.0; | |
| 697 | } | ||
| 698 | |||
| 699 | // --- if new wet event begins then reset IMD & F | ||
| 700 |
2/2✓ Branch 0 taken 1516445 times.
✓ Branch 1 taken 4406 times.
|
1520851 | if ( infil->T <= 0.0 ) |
| 701 | { | ||
| 702 | 1516445 | infil->IMD = (Fumax - infil->Fu) / lu; | |
| 703 | 1516445 | infil->F = 0.0; | |
| 704 | } | ||
| 705 | 1520851 | return 0.0; | |
| 706 | } | ||
| 707 | |||
| 708 | // --- rainfall does not exceed Ksat | ||
| 709 |
2/2✓ Branch 0 taken 805857 times.
✓ Branch 1 taken 7331 times.
|
813188 | if ( ia <= ks ) |
| 710 | { | ||
| 711 | 805857 | dF = ia * tstep; | |
| 712 | 805857 | infil->F += dF; | |
| 713 | 805857 | infil->Fu += dF; | |
| 714 |
2/2✓ Branch 0 taken 437905 times.
✓ Branch 1 taken 367952 times.
|
805857 | infil->Fu = MIN(infil->Fu, Fumax); |
| 715 |
4/4✓ Branch 0 taken 805287 times.
✓ Branch 1 taken 570 times.
✓ Branch 2 taken 711599 times.
✓ Branch 3 taken 93688 times.
|
805857 | if ( modelType == GREEN_AMPT && infil->T <= 0.0 ) |
| 716 | { | ||
| 717 | 711599 | infil->IMD = (Fumax - infil->Fu) / lu; | |
| 718 | 711599 | infil->F = 0.0; | |
| 719 | } | ||
| 720 | 805857 | return ia; | |
| 721 | } | ||
| 722 | |||
| 723 | // --- rainfall exceeds Ksat; renew time to drain upper zone | ||
| 724 | 7331 | infil->T = 5400.0 / lu / Evap.recoveryFactor; | |
| 725 | |||
| 726 | // --- find volume needed to saturate surface layer | ||
| 727 | 7331 | Fs = ks * (infil->S + depth) * infil->IMD / (ia - ks); | |
| 728 | |||
| 729 | // --- surface layer already saturated | ||
| 730 |
2/2✓ Branch 0 taken 310 times.
✓ Branch 1 taken 7021 times.
|
7331 | if ( infil->F > Fs ) |
| 731 | { | ||
| 732 | 310 | infil->Sat = TRUE; | |
| 733 | 310 | return grnampt_getSatInfil(infil, tstep, irate, depth); | |
| 734 | } | ||
| 735 | |||
| 736 | // --- surface layer remains unsaturated | ||
| 737 |
2/2✓ Branch 0 taken 6608 times.
✓ Branch 1 taken 413 times.
|
7021 | if ( infil->F + ia*tstep < Fs ) |
| 738 | { | ||
| 739 | 6608 | dF = ia * tstep; | |
| 740 | 6608 | infil->F += dF; | |
| 741 | 6608 | infil->Fu += dF; | |
| 742 |
2/2✓ Branch 0 taken 3343 times.
✓ Branch 1 taken 3265 times.
|
6608 | infil->Fu = MIN(infil->Fu, Fumax); |
| 743 | 6608 | return ia; | |
| 744 | } | ||
| 745 | |||
| 746 | // --- surface layer becomes saturated during time step; | ||
| 747 | // --- compute portion of tstep when saturated | ||
| 748 | 413 | ts = tstep - (Fs - infil->F) / ia; | |
| 749 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 413 times.
|
413 | if ( ts <= 0.0 ) ts = 0.0; |
| 750 | |||
| 751 | // --- compute new total volume infiltrated | ||
| 752 | 413 | c1 = (infil->S + depth) * infil->IMD; | |
| 753 | 413 | F2 = grnampt_getF2(Fs, c1, ks, ts); | |
| 754 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 407 times.
|
413 | if ( F2 > Fs + ia*ts ) F2 = Fs + ia*ts; |
| 755 | |||
| 756 | // --- compute infiltration rate | ||
| 757 | 413 | dF = F2 - infil->F; | |
| 758 | 413 | infil->F = F2; | |
| 759 | 413 | infil->Fu += dF; | |
| 760 |
2/2✓ Branch 0 taken 23 times.
✓ Branch 1 taken 390 times.
|
413 | infil->Fu = MIN(infil->Fu, Fumax); |
| 761 | 413 | infil->Sat = TRUE; | |
| 762 | 413 | return dF / tstep; | |
| 763 | } | ||
| 764 | |||
| 765 | //============================================================================= | ||
| 766 | |||
| 767 | 61432 | double grnampt_getSatInfil(TGrnAmpt *infil, double tstep, double irate, | |
| 768 | double depth) | ||
| 769 | // | ||
| 770 | // Input: infil = ptr. to Green-Ampt infiltration object | ||
| 771 | // tstep = runoff time step (sec), | ||
| 772 | // irate = net "rainfall" rate to upper zone (ft/sec); | ||
| 773 | // = rainfall + snowmelt + runon, | ||
| 774 | // does not include ponded water (added on below) | ||
| 775 | // depth = depth of ponded water (ft). | ||
| 776 | // Output: returns infiltration rate (ft/sec) | ||
| 777 | // Purpose: computes Green-Ampt infiltration when upper soil zone is | ||
| 778 | // saturated. | ||
| 779 | // | ||
| 780 | { | ||
| 781 | double ia, c1, dF, F2; | ||
| 782 | 61432 | double ks = infil->Ks * InfilFactor; | |
| 783 | 61432 | double lu = infil->Lu * sqrt(InfilFactor); | |
| 784 | |||
| 785 | // --- get available infiltration rate (rainfall + ponded water) | ||
| 786 | 61432 | ia = irate + depth / tstep; | |
| 787 |
2/2✓ Branch 0 taken 2584 times.
✓ Branch 1 taken 58848 times.
|
61432 | if ( ia < ZERO ) return 0.0; |
| 788 | |||
| 789 | // --- re-set new event recovery time | ||
| 790 | 58848 | infil->T = 5400.0 / lu / Evap.recoveryFactor; | |
| 791 | |||
| 792 | // --- solve G-A equation for new cumulative infiltration volume (F2) | ||
| 793 | 58848 | c1 = (infil->S + depth) * infil->IMD; | |
| 794 | 58848 | F2 = grnampt_getF2(infil->F, c1, ks, tstep); | |
| 795 | 58848 | dF = F2 - infil->F; | |
| 796 | |||
| 797 | // --- all available water infiltrates -- set saturated state to false | ||
| 798 |
2/2✓ Branch 0 taken 709 times.
✓ Branch 1 taken 58139 times.
|
58848 | if ( dF > ia * tstep ) |
| 799 | { | ||
| 800 | 709 | dF = ia * tstep; | |
| 801 | 709 | infil->Sat = FALSE; | |
| 802 | } | ||
| 803 | |||
| 804 | // --- update total infiltration and upper zone moisture deficit | ||
| 805 | 58848 | infil->F += dF; | |
| 806 | 58848 | infil->Fu += dF; | |
| 807 |
2/2✓ Branch 0 taken 3172 times.
✓ Branch 1 taken 55676 times.
|
58848 | infil->Fu = MIN(infil->Fu, Fumax); |
| 808 | 58848 | return dF / tstep; | |
| 809 | } | ||
| 810 | |||
| 811 | //============================================================================= | ||
| 812 | |||
| 813 | 59261 | double grnampt_getF2(double f1, double c1, double ks, double ts) | |
| 814 | // | ||
| 815 | // Input: f1 = old infiltration volume (ft) | ||
| 816 | // c1 = head * moisture deficit (ft) | ||
| 817 | // ks = sat. hyd. conductivity (ft/sec) | ||
| 818 | // ts = time step (sec) | ||
| 819 | // Output: returns infiltration volume at end of time step (ft) | ||
| 820 | // Purpose: computes new infiltration volume over a time step | ||
| 821 | // using Green-Ampt formula for saturated upper soil zone. | ||
| 822 | // | ||
| 823 | { | ||
| 824 | int i; | ||
| 825 | 59261 | double f2 = f1; | |
| 826 | double f2min; | ||
| 827 | double df2; | ||
| 828 | double c2; | ||
| 829 | |||
| 830 | // --- find min. infil. volume | ||
| 831 | 59261 | f2min = f1 + ks * ts; | |
| 832 | |||
| 833 | // --- use min. infil. volume for 0 moisture deficit | ||
| 834 |
2/2✓ Branch 0 taken 2049 times.
✓ Branch 1 taken 57212 times.
|
59261 | if ( c1 == 0.0 ) return f2min; |
| 835 | |||
| 836 | // --- use direct form of G-A equation for small time steps | ||
| 837 | // and c1/f1 < 100 | ||
| 838 |
4/4✓ Branch 0 taken 36002 times.
✓ Branch 1 taken 21210 times.
✓ Branch 2 taken 35967 times.
✓ Branch 3 taken 35 times.
|
57212 | if ( ts < 10.0 && f1 > 0.01 * c1 ) |
| 839 | { | ||
| 840 | 35967 | f2 = f1 + ks * (1.0 + c1/f1) * ts; | |
| 841 |
1/2✓ Branch 0 taken 35967 times.
✗ Branch 1 not taken.
|
35967 | return MAX(f2, f2min); |
| 842 | } | ||
| 843 | |||
| 844 | // --- use Newton-Raphson method to solve integrated G-A equation | ||
| 845 | // (convergence limit reduced from that used in previous releases) | ||
| 846 | 21245 | c2 = c1 * log(f1 + c1) - ks * ts; | |
| 847 |
1/2✓ Branch 0 taken 49893 times.
✗ Branch 1 not taken.
|
49893 | for ( i = 1; i <= 20; i++ ) |
| 848 | { | ||
| 849 | 49893 | df2 = (f2 - f1 - c1 * log(f2 + c1) + c2) / (1.0 - c1 / (f2 + c1) ); | |
| 850 |
2/2✓ Branch 0 taken 21245 times.
✓ Branch 1 taken 28648 times.
|
49893 | if ( fabs(df2) < 0.00001 ) |
| 851 | { | ||
| 852 |
1/2✓ Branch 0 taken 21245 times.
✗ Branch 1 not taken.
|
21245 | return MAX(f2, f2min); |
| 853 | } | ||
| 854 | 28648 | f2 -= df2; | |
| 855 | } | ||
| 856 | ✗ | return f2min; | |
| 857 | } | ||
| 858 | |||
| 859 | //============================================================================= | ||
| 860 | |||
| 861 | 9 | int curvenum_setParams(TCurveNum *infil, double p[]) | |
| 862 | // | ||
| 863 | // Input: infil = ptr. to Curve Number infiltration object | ||
| 864 | // p[] = array of parameter values | ||
| 865 | // Output: returns TRUE if parameters are valid, FALSE otherwise | ||
| 866 | // Purpose: assigns Curve Number infiltration parameters to a subcatchment. | ||
| 867 | // | ||
| 868 | { | ||
| 869 | |||
| 870 | // --- convert Curve Number to max. infil. capacity | ||
| 871 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
9 | if ( p[0] < 10.0 ) p[0] = 10.0; |
| 872 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
9 | if ( p[0] > 99.0 ) p[0] = 99.0; |
| 873 | 9 | infil->Smax = (1000.0 / p[0] - 10.0) / 12.0; | |
| 874 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
9 | if ( infil->Smax < 0.0 ) return FALSE; |
| 875 | |||
| 876 | // --- convert drying time (days) to a regeneration const. (1/sec) | ||
| 877 |
1/2✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
|
9 | if ( p[2] > 0.0 ) infil->regen = 1.0 / (p[2] * SECperDAY); |
| 878 | ✗ | else return FALSE; | |
| 879 | |||
| 880 | // --- compute inter-event time from regeneration const. as in Green-Ampt | ||
| 881 | 9 | infil->Tmax = 0.06 / infil->regen; | |
| 882 | |||
| 883 | 9 | return TRUE; | |
| 884 | } | ||
| 885 | |||
| 886 | //============================================================================= | ||
| 887 | |||
| 888 | 9 | void curvenum_initState(TCurveNum *infil) | |
| 889 | // | ||
| 890 | // Input: infil = ptr. to Curve Number infiltration object | ||
| 891 | // Output: none | ||
| 892 | // Purpose: initializes state of Curve Number infiltration for a subcatchment. | ||
| 893 | // | ||
| 894 | { | ||
| 895 | 9 | infil->S = infil->Smax; | |
| 896 | 9 | infil->P = 0.0; | |
| 897 | 9 | infil->F = 0.0; | |
| 898 | 9 | infil->T = 0.0; | |
| 899 | 9 | infil->Se = infil->Smax; | |
| 900 | 9 | infil->f = 0.0; | |
| 901 | 9 | } | |
| 902 | |||
| 903 | ✗ | void curvenum_getState(TCurveNum *infil, double x[]) | |
| 904 | { | ||
| 905 | ✗ | x[0] = infil->S; | |
| 906 | ✗ | x[1] = infil->P; | |
| 907 | ✗ | x[2] = infil->F; | |
| 908 | ✗ | x[3] = infil->T; | |
| 909 | ✗ | x[4] = infil->Se; | |
| 910 | ✗ | x[5] = infil->f; | |
| 911 | ✗ | } | |
| 912 | |||
| 913 | ✗ | void curvenum_setState(TCurveNum *infil, double x[]) | |
| 914 | { | ||
| 915 | ✗ | infil->S = x[0]; | |
| 916 | ✗ | infil->P = x[1]; | |
| 917 | ✗ | infil->F = x[2]; | |
| 918 | ✗ | infil->T = x[3]; | |
| 919 | ✗ | infil->Se = x[4]; | |
| 920 | ✗ | infil->f = x[5]; | |
| 921 | ✗ | } | |
| 922 | |||
| 923 | //============================================================================= | ||
| 924 | |||
| 925 | 927 | double curvenum_getInfil(TCurveNum *infil, double tstep, double irate, | |
| 926 | double depth) | ||
| 927 | // | ||
| 928 | // Input: infil = ptr. to Curve Number infiltration object | ||
| 929 | // tstep = runoff time step (sec), | ||
| 930 | // irate = rainfall rate (ft/sec); | ||
| 931 | // depth = depth of runon + ponded water (ft) | ||
| 932 | // Output: returns infiltration rate (ft/sec) | ||
| 933 | // Purpose: computes infiltration rate using the Curve Number method. | ||
| 934 | // Note: this function treats runon from other subcatchments as part | ||
| 935 | // of the ponded depth and not as an effective rainfall rate. | ||
| 936 | { | ||
| 937 | double F1; // new cumulative infiltration (ft) | ||
| 938 | 927 | double f1 = 0.0; // new infiltration rate (ft/sec) | |
| 939 | 927 | double fa = irate + depth/tstep; // max. available infil. rate (ft/sec) | |
| 940 | |||
| 941 | // --- case where there is rainfall | ||
| 942 |
2/2✓ Branch 0 taken 216 times.
✓ Branch 1 taken 711 times.
|
927 | if ( irate > ZERO ) |
| 943 | { | ||
| 944 | // --- check if new rain event | ||
| 945 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 216 times.
|
216 | if ( infil->T >= infil->Tmax ) |
| 946 | { | ||
| 947 | ✗ | infil->P = 0.0; | |
| 948 | ✗ | infil->F = 0.0; | |
| 949 | ✗ | infil->f = 0.0; | |
| 950 | ✗ | infil->Se = infil->S; | |
| 951 | } | ||
| 952 | 216 | infil->T = 0.0; | |
| 953 | |||
| 954 | // --- update cumulative precip. | ||
| 955 | 216 | infil->P += irate * tstep; | |
| 956 | |||
| 957 | // --- find potential new cumulative infiltration | ||
| 958 | 216 | F1 = infil->P * (1.0 - infil->P / (infil->P + infil->Se)); | |
| 959 | |||
| 960 | // --- compute potential infiltration rate | ||
| 961 | 216 | f1 = (F1 - infil->F) / tstep; | |
| 962 |
2/4✓ Branch 0 taken 216 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 216 times.
|
216 | if ( f1 < 0.0 || infil->S <= 0.0 ) f1 = 0.0; |
| 963 | |||
| 964 | } | ||
| 965 | |||
| 966 | // --- case of no rainfall | ||
| 967 | else | ||
| 968 | { | ||
| 969 | // --- if there is ponded water then use previous infil. rate | ||
| 970 |
3/4✓ Branch 0 taken 289 times.
✓ Branch 1 taken 422 times.
✓ Branch 2 taken 289 times.
✗ Branch 3 not taken.
|
711 | if ( depth > MIN_TOTAL_DEPTH && infil->S > 0.0 ) |
| 971 | { | ||
| 972 | 289 | f1 = infil->f; | |
| 973 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 289 times.
|
289 | if ( f1*tstep > infil->S ) f1 = infil->S / tstep; |
| 974 | } | ||
| 975 | |||
| 976 | // --- otherwise update inter-event time | ||
| 977 | 422 | else infil->T += tstep; | |
| 978 | } | ||
| 979 | |||
| 980 | // --- if there is some infiltration | ||
| 981 |
2/2✓ Branch 0 taken 505 times.
✓ Branch 1 taken 422 times.
|
927 | if ( f1 > 0.0 ) |
| 982 | { | ||
| 983 | // --- limit infil. rate to max. available rate | ||
| 984 |
1/2✓ Branch 0 taken 505 times.
✗ Branch 1 not taken.
|
505 | f1 = MIN(f1, fa); |
| 985 |
1/2✓ Branch 0 taken 505 times.
✗ Branch 1 not taken.
|
505 | f1 = MAX(f1, 0.0); |
| 986 | |||
| 987 | // --- update actual cumulative infiltration | ||
| 988 | 505 | infil->F += f1 * tstep; | |
| 989 | |||
| 990 | // --- reduce infil. capacity if a regen. constant was supplied | ||
| 991 |
1/2✓ Branch 0 taken 505 times.
✗ Branch 1 not taken.
|
505 | if ( infil->regen > 0.0 ) |
| 992 | { | ||
| 993 | 505 | infil->S -= f1 * tstep; | |
| 994 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 505 times.
|
505 | if ( infil->S < 0.0 ) infil->S = 0.0; |
| 995 | } | ||
| 996 | } | ||
| 997 | |||
| 998 | // --- otherwise regenerate infil. capacity | ||
| 999 | else | ||
| 1000 | { | ||
| 1001 | 422 | infil->S += infil->regen * infil->Smax * tstep * Evap.recoveryFactor; | |
| 1002 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 422 times.
|
422 | if ( infil->S > infil->Smax ) infil->S = infil->Smax; |
| 1003 | } | ||
| 1004 | 927 | infil->f = f1; | |
| 1005 | 927 | return f1; | |
| 1006 | } | ||
| 1007 |