gwater.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // gwater.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 | // Groundwater functions. | ||
| 10 | // | ||
| 11 | // Update History | ||
| 12 | // ============== | ||
| 13 | // Build 5.1.007: | ||
| 14 | // - User-supplied function for deep GW seepage flow added. | ||
| 15 | // - New variable names for use in user-supplied GW flow equations added. | ||
| 16 | // Build 5.1.008: | ||
| 17 | // - More variable names for user-supplied GW flow equations added. | ||
| 18 | // - Subcatchment area made into a shared variable. | ||
| 19 | // - Evaporation loss initialized to 0. | ||
| 20 | // - Support for collecting GW statistics added. | ||
| 21 | // Build 5.1.010: | ||
| 22 | // - Unsaturated hydraulic conductivity added to GW flow equation variables. | ||
| 23 | //----------------------------------------------------------------------------- | ||
| 24 | #define _CRT_SECURE_NO_DEPRECATE | ||
| 25 | |||
| 26 | #include <stdlib.h> | ||
| 27 | #include <string.h> | ||
| 28 | #include <math.h> | ||
| 29 | #include "headers.h" | ||
| 30 | #include "odesolve.h" | ||
| 31 | |||
| 32 | //----------------------------------------------------------------------------- | ||
| 33 | // Constants | ||
| 34 | //----------------------------------------------------------------------------- | ||
| 35 | static const double GWTOL = 0.0001; // ODE solver tolerance | ||
| 36 | static const double XTOL = 0.001; // tolerance for moisture & depth | ||
| 37 | |||
| 38 | enum GWstates {THETA, // moisture content of upper GW zone | ||
| 39 | LOWERDEPTH}; // depth of lower saturated GW zone | ||
| 40 | |||
| 41 | enum GWvariables { | ||
| 42 | gwvHGW, // water table height (ft) | ||
| 43 | gwvHSW, // surface water height (ft) | ||
| 44 | gwvHCB, // channel bottom height (ft) | ||
| 45 | gwvHGS, // ground surface height (ft) | ||
| 46 | gwvKS, // sat. hyd. condutivity (ft/s) | ||
| 47 | gwvK, // unsat. hyd. conductivity (ft/s) | ||
| 48 | gwvTHETA, // upper zone moisture content | ||
| 49 | gwvPHI, // soil porosity | ||
| 50 | gwvFI, // surface infiltration (ft/s) | ||
| 51 | gwvFU, // uper zone percolation rate (ft/s) | ||
| 52 | gwvA, // subcatchment area (ft2) | ||
| 53 | gwvMAX}; | ||
| 54 | |||
| 55 | // Names of GW variables that can be used in GW outflow expression | ||
| 56 | static char* GWVarWords[] = {"HGW", "HSW", "HCB", "HGS", "KS", "K", | ||
| 57 | "THETA", "PHI", "FI", "FU", "A", NULL}; | ||
| 58 | |||
| 59 | //----------------------------------------------------------------------------- | ||
| 60 | // Shared variables | ||
| 61 | //----------------------------------------------------------------------------- | ||
| 62 | // NOTE: all flux rates are in ft/sec, all depths are in ft. | ||
| 63 | static double Area; // subcatchment area (ft2) | ||
| 64 | static double Infil; // infiltration rate from surface | ||
| 65 | static double MaxEvap; // max. evaporation rate | ||
| 66 | static double AvailEvap; // available evaporation rate | ||
| 67 | static double UpperEvap; // evaporation rate from upper GW zone | ||
| 68 | static double LowerEvap; // evaporation rate from lower GW zone | ||
| 69 | static double UpperPerc; // percolation rate from upper to lower zone | ||
| 70 | static double LowerLoss; // loss rate from lower GW zone | ||
| 71 | static double GWFlow; // flow rate from lower zone to conveyance node | ||
| 72 | static double MaxUpperPerc; // upper limit on UpperPerc | ||
| 73 | static double MaxGWFlowPos; // upper limit on GWFlow when its positve | ||
| 74 | static double MaxGWFlowNeg; // upper limit on GWFlow when its negative | ||
| 75 | static double FracPerv; // fraction of surface that is pervious | ||
| 76 | static double TotalDepth; // total depth of GW aquifer | ||
| 77 | static double Theta; // moisture content of upper zone | ||
| 78 | static double HydCon; // unsaturated hydraulic conductivity (ft/s) | ||
| 79 | static double Hgw; // ht. of saturated zone | ||
| 80 | static double Hstar; // ht. from aquifer bottom to node invert | ||
| 81 | static double Hsw; // ht. from aquifer bottom to water surface | ||
| 82 | static double Tstep; // current time step (sec) | ||
| 83 | static TAquifer A; // aquifer being analyzed | ||
| 84 | static TGroundwater* GW; // groundwater object being analyzed | ||
| 85 | static MathExpr* LatFlowExpr; // user-supplied lateral GW flow expression | ||
| 86 | static MathExpr* DeepFlowExpr; // user-supplied deep GW flow expression | ||
| 87 | |||
| 88 | //----------------------------------------------------------------------------- | ||
| 89 | // External Functions (declared in funcs.h) | ||
| 90 | //----------------------------------------------------------------------------- | ||
| 91 | // gwater_readAquiferParams (called by input_readLine) | ||
| 92 | // gwater_readGroundwaterParams (called by input_readLine) | ||
| 93 | // gwater_readFlowExpression (called by input_readLine) | ||
| 94 | // gwater_deleteFlowExpression (called by deleteObjects in project.c) | ||
| 95 | // gwater_validateAquifer (called by swmm_open) | ||
| 96 | // gwater_validate (called by subcatch_validate) | ||
| 97 | // gwater_initState (called by subcatch_initState) | ||
| 98 | // gwater_getVolume (called by massbal_open & massbal_getGwaterError) | ||
| 99 | // gwater_getGroundwater (called by getSubareaRunoff in subcatch.c) | ||
| 100 | // gwater_getState (called by saveRunoff in hotstart.c) | ||
| 101 | // gwater_setState (called by readRunoff in hotstart.c) | ||
| 102 | |||
| 103 | //----------------------------------------------------------------------------- | ||
| 104 | // Local functions | ||
| 105 | //----------------------------------------------------------------------------- | ||
| 106 | static void getDxDt(double t, double* x, double* dxdt); | ||
| 107 | static void getFluxes(double upperVolume, double lowerDepth); | ||
| 108 | static void getEvapRates(double theta, double upperDepth); | ||
| 109 | static double getUpperPerc(double theta, double upperDepth); | ||
| 110 | static double getGWFlow(double lowerDepth); | ||
| 111 | static void updateMassBal(double area, double tStep); | ||
| 112 | |||
| 113 | // Used to process custom GW outflow equations | ||
| 114 | static int getVariableIndex(char* s); | ||
| 115 | static double getVariableValue(int varIndex); | ||
| 116 | |||
| 117 | //============================================================================= | ||
| 118 | |||
| 119 | 2186 | int gwater_readAquiferParams(int j, char* tok[], int ntoks) | |
| 120 | // | ||
| 121 | // Input: j = aquifer index | ||
| 122 | // tok[] = array of string tokens | ||
| 123 | // ntoks = number of tokens | ||
| 124 | // Output: returns error message | ||
| 125 | // Purpose: reads aquifer parameter values from line of input data | ||
| 126 | // | ||
| 127 | // Data line contains following parameters: | ||
| 128 | // ID, porosity, wiltingPoint, fieldCapacity, conductivity, | ||
| 129 | // conductSlope, tensionSlope, upperEvapFraction, lowerEvapDepth, | ||
| 130 | // gwRecession, bottomElev, waterTableElev, upperMoisture | ||
| 131 | // (evapPattern) | ||
| 132 | // | ||
| 133 | { | ||
| 134 | int i, p; | ||
| 135 | double x[12]; | ||
| 136 | char *id; | ||
| 137 | |||
| 138 | // --- check that aquifer exists | ||
| 139 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2186 times.
|
2186 | if ( ntoks < 13 ) return error_setInpError(ERR_ITEMS, ""); |
| 140 | 2186 | id = project_findID(AQUIFER, tok[0]); | |
| 141 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2186 times.
|
2186 | if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]); |
| 142 | |||
| 143 | // --- read remaining tokens as numbers | ||
| 144 |
2/2✓ Branch 0 taken 24046 times.
✓ Branch 1 taken 2186 times.
|
26232 | for (i = 0; i < 11; i++) x[i] = 0.0; |
| 145 |
2/2✓ Branch 0 taken 26232 times.
✓ Branch 1 taken 2186 times.
|
28418 | for (i = 1; i < 13; i++) |
| 146 | { | ||
| 147 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 26232 times.
|
26232 | if ( ! getDouble(tok[i], &x[i-1]) ) |
| 148 | ✗ | return error_setInpError(ERR_NUMBER, tok[i]); | |
| 149 | } | ||
| 150 | |||
| 151 | // --- read upper evap pattern if present | ||
| 152 | 2186 | p = -1; | |
| 153 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2186 times.
|
2186 | if ( ntoks > 13 ) |
| 154 | { | ||
| 155 | ✗ | p = project_findObject(TIMEPATTERN, tok[13]); | |
| 156 | ✗ | if ( p < 0 ) return error_setInpError(ERR_NAME, tok[13]); | |
| 157 | } | ||
| 158 | |||
| 159 | // --- assign parameters to aquifer object | ||
| 160 | 2186 | Aquifer[j].ID = id; | |
| 161 | 2186 | Aquifer[j].porosity = x[0]; | |
| 162 | 2186 | Aquifer[j].wiltingPoint = x[1]; | |
| 163 | 2186 | Aquifer[j].fieldCapacity = x[2]; | |
| 164 | 2186 | Aquifer[j].conductivity = x[3] / UCF(RAINFALL); | |
| 165 | 2186 | Aquifer[j].conductSlope = x[4]; | |
| 166 | 2186 | Aquifer[j].tensionSlope = x[5] / UCF(LENGTH); | |
| 167 | 2186 | Aquifer[j].upperEvapFrac = x[6]; | |
| 168 | 2186 | Aquifer[j].lowerEvapDepth = x[7] / UCF(LENGTH); | |
| 169 | 2186 | Aquifer[j].lowerLossCoeff = x[8] / UCF(RAINFALL); | |
| 170 | 2186 | Aquifer[j].bottomElev = x[9] / UCF(LENGTH); | |
| 171 | 2186 | Aquifer[j].waterTableElev = x[10] / UCF(LENGTH); | |
| 172 | 2186 | Aquifer[j].upperMoisture = x[11]; | |
| 173 | 2186 | Aquifer[j].upperEvapPat = p; | |
| 174 | 2186 | return 0; | |
| 175 | } | ||
| 176 | |||
| 177 | //============================================================================= | ||
| 178 | |||
| 179 | 2241 | int gwater_readGroundwaterParams(char* tok[], int ntoks) | |
| 180 | // | ||
| 181 | // Input: tok[] = array of string tokens | ||
| 182 | // ntoks = number of tokens | ||
| 183 | // Output: returns error code | ||
| 184 | // Purpose: reads groundwater inflow parameters for a subcatchment from | ||
| 185 | // a line of input data. | ||
| 186 | // | ||
| 187 | // Data format is: | ||
| 188 | // subcatch aquifer node surfElev a1 b1 a2 b2 a3 fixedDepth + | ||
| 189 | // (nodeElev bottomElev waterTableElev upperMoisture ) | ||
| 190 | // | ||
| 191 | { | ||
| 192 | int i, j, k, m, n; | ||
| 193 | double x[11]; | ||
| 194 | TGroundwater* gw; | ||
| 195 | |||
| 196 | // --- check that specified subcatchment, aquifer & node exist | ||
| 197 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, ""); |
| 198 | 2241 | j = project_findObject(SUBCATCH, tok[0]); | |
| 199 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]); |
| 200 | |||
| 201 | // --- check for enough tokens | ||
| 202 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( ntoks < 11 ) return error_setInpError(ERR_ITEMS, ""); |
| 203 | |||
| 204 | // --- check that specified aquifer and node exists | ||
| 205 | 2241 | k = project_findObject(AQUIFER, tok[1]); | |
| 206 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( k < 0 ) return error_setInpError(ERR_NAME, tok[1]); |
| 207 | 2241 | n = project_findObject(NODE, tok[2]); | |
| 208 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( n < 0 ) return error_setInpError(ERR_NAME, tok[2]); |
| 209 | |||
| 210 | // -- read in the flow parameters | ||
| 211 |
2/2✓ Branch 0 taken 15687 times.
✓ Branch 1 taken 2241 times.
|
17928 | for ( i = 0; i < 7; i++ ) |
| 212 | { | ||
| 213 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 15687 times.
|
15687 | if ( ! getDouble(tok[i+3], &x[i]) ) |
| 214 | ✗ | return error_setInpError(ERR_NUMBER, tok[i+3]); | |
| 215 | } | ||
| 216 | |||
| 217 | // --- read in optional depth parameters | ||
| 218 |
2/2✓ Branch 0 taken 8964 times.
✓ Branch 1 taken 2241 times.
|
11205 | for ( i = 7; i < 11; i++) |
| 219 | { | ||
| 220 | 8964 | x[i] = MISSING; | |
| 221 | 8964 | m = i + 3; | |
| 222 |
3/4✓ Branch 0 taken 8964 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 6731 times.
✓ Branch 3 taken 2233 times.
|
8964 | if ( ntoks > m && *tok[m] != '*' ) |
| 223 | { | ||
| 224 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 6731 times.
|
6731 | if (! getDouble(tok[m], &x[i]) ) |
| 225 | ✗ | return error_setInpError(ERR_NUMBER, tok[m]); | |
| 226 |
2/2✓ Branch 0 taken 4490 times.
✓ Branch 1 taken 2241 times.
|
6731 | if ( i < 10 ) x[i] /= UCF(LENGTH); |
| 227 | } | ||
| 228 | } | ||
| 229 | |||
| 230 | // --- create a groundwater flow object | ||
| 231 |
1/2✓ Branch 0 taken 2241 times.
✗ Branch 1 not taken.
|
2241 | if ( !Subcatch[j].groundwater ) |
| 232 | { | ||
| 233 | 2241 | gw = (TGroundwater *) malloc(sizeof(TGroundwater)); | |
| 234 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( !gw ) return error_setInpError(ERR_MEMORY, ""); |
| 235 | 2241 | Subcatch[j].groundwater = gw; | |
| 236 | } | ||
| 237 | ✗ | else gw = Subcatch[j].groundwater; | |
| 238 | |||
| 239 | // --- populate the groundwater flow object with its parameters | ||
| 240 | 2241 | gw->aquifer = k; | |
| 241 | 2241 | gw->node = n; | |
| 242 | 2241 | gw->surfElev = x[0] / UCF(LENGTH); | |
| 243 | 2241 | gw->a1 = x[1]; | |
| 244 | 2241 | gw->b1 = x[2]; | |
| 245 | 2241 | gw->a2 = x[3]; | |
| 246 | 2241 | gw->b2 = x[4]; | |
| 247 | 2241 | gw->a3 = x[5]; | |
| 248 | 2241 | gw->fixedDepth = x[6] / UCF(LENGTH); | |
| 249 | 2241 | gw->nodeElev = x[7]; //already converted to ft. | |
| 250 | 2241 | gw->bottomElev = x[8]; | |
| 251 | 2241 | gw->waterTableElev = x[9]; | |
| 252 | 2241 | gw->upperMoisture = x[10]; | |
| 253 | 2241 | return 0; | |
| 254 | } | ||
| 255 | |||
| 256 | //============================================================================= | ||
| 257 | |||
| 258 | 8 | int gwater_readFlowExpression(char* tok[], int ntoks) | |
| 259 | // | ||
| 260 | // Input: tok[] = array of string tokens | ||
| 261 | // ntoks = number of tokens | ||
| 262 | // Output: returns error code | ||
| 263 | // Purpose: reads mathematical expression for lateral or deep groundwater | ||
| 264 | // flow for a subcatchment from a line of input data. | ||
| 265 | // | ||
| 266 | // Format is: subcatch LATERAL/DEEP <expr> | ||
| 267 | // where subcatch is the ID of the subcatchment, LATERAL is for lateral | ||
| 268 | // GW flow, DEEP is for deep GW flow and <expr> is any well-formed math | ||
| 269 | // expression. | ||
| 270 | // | ||
| 271 | { | ||
| 272 | int i, j, k; | ||
| 273 | char exprStr[MAXLINE+1]; | ||
| 274 | MathExpr* expr; | ||
| 275 | |||
| 276 | // --- return if too few tokens | ||
| 277 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, ""); |
| 278 | |||
| 279 | // --- check that subcatchment exists | ||
| 280 | 8 | j = project_findObject(SUBCATCH, tok[0]); | |
| 281 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]); |
| 282 | |||
| 283 | // --- check if expression is for lateral or deep GW flow | ||
| 284 | 8 | k = 1; | |
| 285 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
|
8 | if ( match(tok[1], "LAT") ) k = 1; |
| 286 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | else if ( match(tok[1], "DEEP") ) k = 2; |
| 287 | ✗ | else return error_setInpError(ERR_KEYWORD, tok[1]); | |
| 288 | |||
| 289 | // --- concatenate remaining tokens into a single string | ||
| 290 | 8 | sstrncpy(exprStr, tok[2], MAXLINE); | |
| 291 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | for ( i = 3; i < ntoks; i++) |
| 292 | { | ||
| 293 | ✗ | sstrcat(exprStr, " ", MAXLINE+1); | |
| 294 | ✗ | sstrcat(exprStr, tok[i], MAXLINE+1); | |
| 295 | } | ||
| 296 | |||
| 297 | // --- delete any previous flow eqn. | ||
| 298 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if ( k == 1 ) mathexpr_delete(Subcatch[j].gwLatFlowExpr); |
| 299 | 8 | else mathexpr_delete(Subcatch[j].gwDeepFlowExpr); | |
| 300 | |||
| 301 | // --- create a parsed expression tree from the string expr | ||
| 302 | // (getVariableIndex is the function that converts a GW | ||
| 303 | // variable's name into an index number) | ||
| 304 | 8 | expr = mathexpr_create(exprStr, getVariableIndex); | |
| 305 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if ( expr == NULL ) return error_setInpError(ERR_MATH_EXPR, ""); |
| 306 | |||
| 307 | // --- save expression tree with the subcatchment | ||
| 308 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if ( k == 1 ) Subcatch[j].gwLatFlowExpr = expr; |
| 309 | 8 | else Subcatch[j].gwDeepFlowExpr = expr; | |
| 310 | 8 | return 0; | |
| 311 | } | ||
| 312 | |||
| 313 | //============================================================================= | ||
| 314 | |||
| 315 | 2393 | void gwater_deleteFlowExpression(int j) | |
| 316 | // | ||
| 317 | // Input: j = subcatchment index | ||
| 318 | // Output: none | ||
| 319 | // Purpose: deletes a subcatchment's custom groundwater flow expressions. | ||
| 320 | // | ||
| 321 | { | ||
| 322 | 2393 | mathexpr_delete(Subcatch[j].gwLatFlowExpr); | |
| 323 | 2393 | mathexpr_delete(Subcatch[j].gwDeepFlowExpr); | |
| 324 | 2393 | } | |
| 325 | |||
| 326 | //============================================================================= | ||
| 327 | |||
| 328 | 2186 | void gwater_validateAquifer(int j) | |
| 329 | // | ||
| 330 | // Input: j = aquifer index | ||
| 331 | // Output: none | ||
| 332 | // Purpose: validates groundwater aquifer properties . | ||
| 333 | // | ||
| 334 | { | ||
| 335 | int p; | ||
| 336 | |||
| 337 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | if ( Aquifer[j].porosity <= 0.0 |
| 338 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].fieldCapacity >= Aquifer[j].porosity |
| 339 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].wiltingPoint >= Aquifer[j].fieldCapacity |
| 340 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].conductivity <= 0.0 |
| 341 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].conductSlope < 0.0 |
| 342 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].tensionSlope < 0.0 |
| 343 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].upperEvapFrac < 0.0 |
| 344 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].lowerEvapDepth < 0.0 |
| 345 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].waterTableElev < Aquifer[j].bottomElev |
| 346 |
1/2✓ Branch 0 taken 2186 times.
✗ Branch 1 not taken.
|
2186 | || Aquifer[j].upperMoisture > Aquifer[j].porosity |
| 347 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2186 times.
|
2186 | || Aquifer[j].upperMoisture < Aquifer[j].wiltingPoint ) |
| 348 | ✗ | report_writeErrorMsg(ERR_AQUIFER_PARAMS, Aquifer[j].ID); | |
| 349 | |||
| 350 | 2186 | p = Aquifer[j].upperEvapPat; | |
| 351 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 2186 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2186 | if ( p >= 0 && Pattern[p].type != MONTHLY_PATTERN ) |
| 352 | { | ||
| 353 | ✗ | report_writeErrorMsg(ERR_AQUIFER_PARAMS, Aquifer[j].ID); | |
| 354 | } | ||
| 355 | 2186 | } | |
| 356 | |||
| 357 | //============================================================================= | ||
| 358 | |||
| 359 | 2393 | void gwater_validate(int j) | |
| 360 | { | ||
| 361 | TAquifer a; // Aquifer data structure | ||
| 362 | TGroundwater* gw; // Groundwater data structure | ||
| 363 | |||
| 364 | 2393 | gw = Subcatch[j].groundwater; | |
| 365 |
2/2✓ Branch 0 taken 2241 times.
✓ Branch 1 taken 152 times.
|
2393 | if ( gw ) |
| 366 | { | ||
| 367 | 2241 | a = Aquifer[gw->aquifer]; | |
| 368 | |||
| 369 | // ... use aquifer values for missing groundwater parameters | ||
| 370 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( gw->bottomElev == MISSING ) gw->bottomElev = a.bottomElev; |
| 371 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( gw->waterTableElev == MISSING ) gw->waterTableElev = a.waterTableElev; |
| 372 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( gw->upperMoisture == MISSING ) gw->upperMoisture = a.upperMoisture; |
| 373 | |||
| 374 | // ... ground elevation can't be below water table elevation | ||
| 375 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( gw->surfElev < gw->waterTableElev ) |
| 376 | ✗ | report_writeErrorMsg(ERR_GROUND_ELEV, Subcatch[j].ID); | |
| 377 | } | ||
| 378 | 2393 | } | |
| 379 | |||
| 380 | //============================================================================= | ||
| 381 | |||
| 382 | 2241 | void gwater_initState(int j) | |
| 383 | // | ||
| 384 | // Input: j = subcatchment index | ||
| 385 | // Output: none | ||
| 386 | // Purpose: initializes state of subcatchment's groundwater. | ||
| 387 | // | ||
| 388 | { | ||
| 389 | TAquifer a; // Aquifer data structure | ||
| 390 | TGroundwater* gw; // Groundwater data structure | ||
| 391 | |||
| 392 | 2241 | gw = Subcatch[j].groundwater; | |
| 393 |
1/2✓ Branch 0 taken 2241 times.
✗ Branch 1 not taken.
|
2241 | if ( gw ) |
| 394 | { | ||
| 395 | 2241 | a = Aquifer[gw->aquifer]; | |
| 396 | |||
| 397 | // ... initial moisture content | ||
| 398 | 2241 | gw->theta = gw->upperMoisture; | |
| 399 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2233 times.
|
2241 | if ( gw->theta >= a.porosity ) |
| 400 | { | ||
| 401 | 8 | gw->theta = a.porosity - XTOL; | |
| 402 | } | ||
| 403 | |||
| 404 | // ... initial depth of lower (saturated) zone | ||
| 405 | 2241 | gw->lowerDepth = gw->waterTableElev - gw->bottomElev; | |
| 406 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2241 times.
|
2241 | if ( gw->lowerDepth >= gw->surfElev - gw->bottomElev ) |
| 407 | { | ||
| 408 | ✗ | gw->lowerDepth = gw->surfElev - gw->bottomElev - XTOL; | |
| 409 | } | ||
| 410 | |||
| 411 | // ... initial lateral groundwater outflow | ||
| 412 | 2241 | gw->oldFlow = 0.0; | |
| 413 | 2241 | gw->newFlow = 0.0; | |
| 414 | 2241 | gw->evapLoss = 0.0; | |
| 415 | |||
| 416 | // ... initial available infiltration volume into upper zone | ||
| 417 | 4482 | gw->maxInfilVol = (gw->surfElev - gw->waterTableElev) * | |
| 418 | 2241 | (a.porosity - gw->theta) / | |
| 419 | 2241 | subcatch_getFracPerv(j); | |
| 420 | } | ||
| 421 | 2241 | } | |
| 422 | |||
| 423 | //============================================================================= | ||
| 424 | |||
| 425 | 319 | void gwater_getState(int j, double x[]) | |
| 426 | // | ||
| 427 | // Input: j = subcatchment index | ||
| 428 | // Output: x[] = array of groundwater state variables | ||
| 429 | // Purpose: retrieves state of subcatchment's groundwater. | ||
| 430 | // | ||
| 431 | { | ||
| 432 | 319 | TGroundwater* gw = Subcatch[j].groundwater; | |
| 433 | 319 | x[0] = gw->theta; | |
| 434 | 319 | x[1] = gw->bottomElev + gw->lowerDepth; | |
| 435 | 319 | x[2] = gw->newFlow; | |
| 436 | 319 | x[3] = gw->maxInfilVol; | |
| 437 | 319 | } | |
| 438 | |||
| 439 | //============================================================================= | ||
| 440 | |||
| 441 | 319 | void gwater_setState(int j, double x[]) | |
| 442 | // | ||
| 443 | // Input: j = subcatchment index | ||
| 444 | // x[] = array of groundwater state variables | ||
| 445 | // Purpose: assigns values to a subcatchment's groundwater state. | ||
| 446 | // | ||
| 447 | { | ||
| 448 | 319 | TGroundwater* gw = Subcatch[j].groundwater; | |
| 449 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 319 times.
|
319 | if ( gw == NULL ) return; |
| 450 | 319 | gw->theta = x[0]; | |
| 451 | 319 | gw->lowerDepth = x[1] - gw->bottomElev; | |
| 452 | 319 | gw->oldFlow = x[2]; | |
| 453 |
1/2✓ Branch 0 taken 319 times.
✗ Branch 1 not taken.
|
319 | if ( x[3] != MISSING ) gw->maxInfilVol = x[3]; |
| 454 | } | ||
| 455 | |||
| 456 | //============================================================================= | ||
| 457 | |||
| 458 | 4704 | double gwater_getVolume(int j) | |
| 459 | // | ||
| 460 | // Input: j = subcatchment index | ||
| 461 | // Output: returns total volume of groundwater in ft/ft2 | ||
| 462 | // Purpose: finds volume of groundwater stored in upper & lower zones | ||
| 463 | // | ||
| 464 | { | ||
| 465 | TAquifer a; | ||
| 466 | TGroundwater* gw; | ||
| 467 | double upperDepth; | ||
| 468 | 4704 | gw = Subcatch[j].groundwater; | |
| 469 |
2/2✓ Branch 0 taken 222 times.
✓ Branch 1 taken 4482 times.
|
4704 | if ( gw == NULL ) return 0.0; |
| 470 | 4482 | a = Aquifer[gw->aquifer]; | |
| 471 | 4482 | upperDepth = gw->surfElev - gw->bottomElev - gw->lowerDepth; | |
| 472 | 4482 | return (upperDepth * gw->theta) + (gw->lowerDepth * a.porosity); | |
| 473 | } | ||
| 474 | |||
| 475 | //============================================================================= | ||
| 476 | |||
| 477 | 3043013 | void gwater_getGroundwater(int j, double evap, double infil, double tStep) | |
| 478 | // | ||
| 479 | // Purpose: computes groundwater flow from subcatchment during current time step. | ||
| 480 | // Input: j = subcatchment index | ||
| 481 | // evap = pervious surface evaporation volume consumed (ft3) | ||
| 482 | // infil = surface infiltration volume (ft3) | ||
| 483 | // tStep = time step (sec) | ||
| 484 | // Output: none | ||
| 485 | // | ||
| 486 | { | ||
| 487 | int n; // node exchanging groundwater | ||
| 488 | double x[2]; // upper moisture content & lower depth | ||
| 489 | double vUpper; // upper vol. available for percolation | ||
| 490 | double nodeFlow; // max. possible GW flow from node | ||
| 491 | |||
| 492 | // --- save subcatchment's groundwater and aquifer objects to | ||
| 493 | // shared variables | ||
| 494 | 3043013 | GW = Subcatch[j].groundwater; | |
| 495 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3043013 times.
|
3052480 | if ( GW == NULL ) return; |
| 496 | 3043013 | LatFlowExpr = Subcatch[j].gwLatFlowExpr; | |
| 497 | 3043013 | DeepFlowExpr = Subcatch[j].gwDeepFlowExpr; | |
| 498 | 3043013 | A = Aquifer[GW->aquifer]; | |
| 499 | |||
| 500 | // --- get fraction of total area that is pervious | ||
| 501 | 3043013 | FracPerv = subcatch_getFracPerv(j); | |
| 502 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3043013 times.
|
3043013 | if ( FracPerv <= 0.0 ) return; |
| 503 | 3043013 | Area = Subcatch[j].area; | |
| 504 | |||
| 505 | // --- convert infiltration volume (ft3) to equivalent rate | ||
| 506 | // over entire GW (subcatchment) area | ||
| 507 | 3043013 | infil = infil / Area / tStep; | |
| 508 | 3043013 | Infil = infil; | |
| 509 | 3043013 | Tstep = tStep; | |
| 510 | |||
| 511 | // --- convert pervious surface evaporation already exerted (ft3) | ||
| 512 | // to equivalent rate over entire GW (subcatchment) area | ||
| 513 | 3043013 | evap = evap / Area / tStep; | |
| 514 | |||
| 515 | // --- convert max. surface evap rate (ft/sec) to a rate | ||
| 516 | // that applies to GW evap (GW evap can only occur | ||
| 517 | // through the pervious land surface area) | ||
| 518 | 3043013 | MaxEvap = Evap.rate * FracPerv; | |
| 519 | |||
| 520 | // --- available subsurface evaporation is difference between max. | ||
| 521 | // rate and pervious surface evap already exerted | ||
| 522 |
2/2✓ Branch 0 taken 1888070 times.
✓ Branch 1 taken 1154943 times.
|
3043013 | AvailEvap = MAX((MaxEvap - evap), 0.0); |
| 523 | |||
| 524 | // --- save total depth & outlet node properties to shared variables | ||
| 525 | 3043013 | TotalDepth = GW->surfElev - GW->bottomElev; | |
| 526 |
2/2✓ Branch 0 taken 9467 times.
✓ Branch 1 taken 3033546 times.
|
3043013 | if ( TotalDepth <= 0.0 ) return; |
| 527 | 3033546 | n = GW->node; | |
| 528 | |||
| 529 | // --- establish min. water table height above aquifer bottom at which | ||
| 530 | // GW flow can occur (override node's invert if a value was provided | ||
| 531 | // in the GW object) | ||
| 532 |
2/2✓ Branch 0 taken 23040 times.
✓ Branch 1 taken 3010506 times.
|
3033546 | if ( GW->nodeElev != MISSING ) Hstar = GW->nodeElev - GW->bottomElev; |
| 533 | 3010506 | else Hstar = Node[n].invertElev - GW->bottomElev; | |
| 534 | |||
| 535 | // --- establish surface water height (relative to aquifer bottom) | ||
| 536 | // for drainage system node connected to the GW aquifer | ||
| 537 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3033546 times.
|
3033546 | if ( GW->fixedDepth > 0.0 ) |
| 538 | { | ||
| 539 | ✗ | Hsw = GW->fixedDepth + Node[n].invertElev - GW->bottomElev; | |
| 540 | } | ||
| 541 | 3033546 | else Hsw = Node[n].newDepth + Node[n].invertElev - GW->bottomElev; | |
| 542 | |||
| 543 | // --- store state variables (upper zone moisture content, lower zone | ||
| 544 | // depth) in work vector x | ||
| 545 | 3033546 | x[THETA] = GW->theta; | |
| 546 | 3033546 | x[LOWERDEPTH] = GW->lowerDepth; | |
| 547 | |||
| 548 | // --- set limit on percolation rate from upper to lower GW zone | ||
| 549 | 3033546 | vUpper = (TotalDepth - x[LOWERDEPTH]) * (x[THETA] - A.fieldCapacity); | |
| 550 |
2/2✓ Branch 0 taken 2405045 times.
✓ Branch 1 taken 628501 times.
|
3033546 | vUpper = MAX(0.0, vUpper); |
| 551 | 3033546 | MaxUpperPerc = vUpper / tStep; | |
| 552 | |||
| 553 | // --- set limit on GW flow out of aquifer based on volume of lower zone | ||
| 554 | 3033546 | MaxGWFlowPos = x[LOWERDEPTH]*A.porosity / tStep; | |
| 555 | |||
| 556 | // --- set limit on GW flow into aquifer from drainage system node | ||
| 557 | // based on min. of capacity of upper zone and drainage system | ||
| 558 | // inflow to the node | ||
| 559 | 3033546 | MaxGWFlowNeg = (TotalDepth - x[LOWERDEPTH]) * (A.porosity - x[THETA]) | |
| 560 | 3033546 | / tStep; | |
| 561 | 3033546 | nodeFlow = (Node[n].inflow + Node[n].newVolume/tStep) / Area; | |
| 562 |
2/2✓ Branch 0 taken 4674 times.
✓ Branch 1 taken 3028872 times.
|
3033546 | MaxGWFlowNeg = -MIN(MaxGWFlowNeg, nodeFlow); |
| 563 | |||
| 564 | // --- integrate eqns. for d(Theta)/dt and d(LowerDepth)/dt | ||
| 565 | // NOTE: ODE solver must have been initialized previously | ||
| 566 | 3033546 | odesolve_integrate(x, 2, 0, tStep, GWTOL, tStep, getDxDt); | |
| 567 | |||
| 568 | // --- keep state variables within allowable bounds | ||
| 569 |
2/2✓ Branch 0 taken 3033544 times.
✓ Branch 1 taken 2 times.
|
3033546 | x[THETA] = MAX(x[THETA], A.wiltingPoint); |
| 570 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3033546 times.
|
3033546 | if ( x[THETA] >= A.porosity ) |
| 571 | { | ||
| 572 | ✗ | x[THETA] = A.porosity - XTOL; | |
| 573 | ✗ | x[LOWERDEPTH] = TotalDepth - XTOL; | |
| 574 | } | ||
| 575 |
1/2✓ Branch 0 taken 3033546 times.
✗ Branch 1 not taken.
|
3033546 | x[LOWERDEPTH] = MAX(x[LOWERDEPTH], 0.0); |
| 576 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3033546 times.
|
3033546 | if ( x[LOWERDEPTH] >= TotalDepth ) |
| 577 | { | ||
| 578 | ✗ | x[LOWERDEPTH] = TotalDepth - XTOL; | |
| 579 | } | ||
| 580 | |||
| 581 | // --- save new values of state values | ||
| 582 | 3033546 | GW->theta = x[THETA]; | |
| 583 | 3033546 | GW->lowerDepth = x[LOWERDEPTH]; | |
| 584 | 3033546 | getFluxes(GW->theta, GW->lowerDepth); | |
| 585 | 3033546 | GW->oldFlow = GW->newFlow; | |
| 586 | 3033546 | GW->newFlow = GWFlow; | |
| 587 | 3033546 | GW->evapLoss = UpperEvap + LowerEvap; | |
| 588 | |||
| 589 | //--- find max. infiltration volume (as depth over | ||
| 590 | // the pervious portion of the subcatchment) | ||
| 591 | // that upper zone can support in next time step | ||
| 592 | 3033546 | GW->maxInfilVol = (TotalDepth - x[LOWERDEPTH]) * | |
| 593 | 3033546 | (A.porosity - x[THETA]) / FracPerv; | |
| 594 | |||
| 595 | // --- update GW mass balance | ||
| 596 | 3033546 | updateMassBal(Area, tStep); | |
| 597 | |||
| 598 | // --- update GW statistics | ||
| 599 | 3033546 | stats_updateGwaterStats(j, infil, GW->evapLoss, GWFlow, LowerLoss, | |
| 600 | 3033546 | GW->theta, GW->lowerDepth + GW->bottomElev, tStep); | |
| 601 | } | ||
| 602 | |||
| 603 | //============================================================================= | ||
| 604 | |||
| 605 | 3033546 | void updateMassBal(double area, double tStep) | |
| 606 | // | ||
| 607 | // Input: area = subcatchment area (ft2) | ||
| 608 | // tStep = time step (sec) | ||
| 609 | // Output: none | ||
| 610 | // Purpose: updates GW mass balance with volumes of water fluxes. | ||
| 611 | // | ||
| 612 | { | ||
| 613 | double vInfil; // infiltration volume | ||
| 614 | double vUpperEvap; // upper zone evap. volume | ||
| 615 | double vLowerEvap; // lower zone evap. volume | ||
| 616 | double vLowerPerc; // lower zone deep perc. volume | ||
| 617 | double vGwater; // volume of exchanged groundwater | ||
| 618 | 3033546 | double ft2sec = area * tStep; | |
| 619 | |||
| 620 | 3033546 | vInfil = Infil * ft2sec; | |
| 621 | 3033546 | vUpperEvap = UpperEvap * ft2sec; | |
| 622 | 3033546 | vLowerEvap = LowerEvap * ft2sec; | |
| 623 | 3033546 | vLowerPerc = LowerLoss * ft2sec; | |
| 624 | 3033546 | vGwater = 0.5 * (GW->oldFlow + GW->newFlow) * ft2sec; | |
| 625 | 3033546 | massbal_updateGwaterTotals(vInfil, vUpperEvap, vLowerEvap, vLowerPerc, | |
| 626 | vGwater); | ||
| 627 | 3033546 | } | |
| 628 | |||
| 629 | //============================================================================= | ||
| 630 | |||
| 631 | 21234942 | void getFluxes(double theta, double lowerDepth) | |
| 632 | // | ||
| 633 | // Input: upperVolume = vol. depth of upper zone (ft) | ||
| 634 | // upperDepth = depth of upper zone (ft) | ||
| 635 | // Output: none | ||
| 636 | // Purpose: computes water fluxes into/out of upper/lower GW zones. | ||
| 637 | // | ||
| 638 | { | ||
| 639 | double upperDepth; | ||
| 640 | |||
| 641 | // --- find upper zone depth | ||
| 642 |
1/2✓ Branch 0 taken 21234942 times.
✗ Branch 1 not taken.
|
21234942 | lowerDepth = MAX(lowerDepth, 0.0); |
| 643 |
1/2✓ Branch 0 taken 21234942 times.
✗ Branch 1 not taken.
|
21234942 | lowerDepth = MIN(lowerDepth, TotalDepth); |
| 644 | 21234942 | upperDepth = TotalDepth - lowerDepth; | |
| 645 | |||
| 646 | // --- save lower depth and theta to global variables | ||
| 647 | 21234942 | Hgw = lowerDepth; | |
| 648 | 21234942 | Theta = theta; | |
| 649 | |||
| 650 | // --- find evaporation rate from both zones | ||
| 651 | 21234942 | getEvapRates(theta, upperDepth); | |
| 652 | |||
| 653 | // --- find percolation rate from upper to lower zone | ||
| 654 | 21234942 | UpperPerc = getUpperPerc(theta, upperDepth); | |
| 655 |
2/2✓ Branch 0 taken 16990463 times.
✓ Branch 1 taken 4244479 times.
|
21234942 | UpperPerc = MIN(UpperPerc, MaxUpperPerc); |
| 656 | |||
| 657 | // --- find loss rate to deep GW | ||
| 658 |
2/2✓ Branch 0 taken 161400 times.
✓ Branch 1 taken 21073542 times.
|
21234942 | if ( DeepFlowExpr != NULL ) |
| 659 | 322800 | LowerLoss = mathexpr_eval(DeepFlowExpr, getVariableValue) / | |
| 660 | 161400 | UCF(RAINFALL); | |
| 661 | else | ||
| 662 | 21073542 | LowerLoss = A.lowerLossCoeff * lowerDepth / TotalDepth; | |
| 663 |
1/2✓ Branch 0 taken 21234942 times.
✗ Branch 1 not taken.
|
21234942 | LowerLoss = MIN(LowerLoss, lowerDepth/Tstep); |
| 664 | |||
| 665 | // --- find GW flow rate from lower zone to drainage system node | ||
| 666 | 21234942 | GWFlow = getGWFlow(lowerDepth); | |
| 667 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 21234942 times.
|
21234942 | if ( LatFlowExpr != NULL ) |
| 668 | { | ||
| 669 | ✗ | GWFlow += mathexpr_eval(LatFlowExpr, getVariableValue) / UCF(GWFLOW); | |
| 670 | } | ||
| 671 |
3/4✓ Branch 0 taken 21123630 times.
✓ Branch 1 taken 111312 times.
✓ Branch 2 taken 21123630 times.
✗ Branch 3 not taken.
|
21234942 | if ( GWFlow >= 0.0 ) GWFlow = MIN(GWFlow, MaxGWFlowPos); |
| 672 |
2/2✓ Branch 0 taken 5294 times.
✓ Branch 1 taken 106018 times.
|
111312 | else GWFlow = MAX(GWFlow, MaxGWFlowNeg); |
| 673 | 21234942 | } | |
| 674 | |||
| 675 | //============================================================================= | ||
| 676 | |||
| 677 | 18201396 | void getDxDt(double t, double* x, double* dxdt) | |
| 678 | // | ||
| 679 | // Input: t = current time (not used) | ||
| 680 | // x = array of state variables | ||
| 681 | // Output: dxdt = array of time derivatives of state variables | ||
| 682 | // Purpose: computes time derivatives of upper moisture content | ||
| 683 | // and lower depth. | ||
| 684 | // | ||
| 685 | { | ||
| 686 | double qUpper; // inflow - outflow for upper zone (ft/sec) | ||
| 687 | double qLower; // inflow - outflow for lower zone (ft/sec) | ||
| 688 | double denom; | ||
| 689 | |||
| 690 | 18201396 | getFluxes(x[THETA], x[LOWERDEPTH]); | |
| 691 | 18201396 | qUpper = Infil - UpperEvap - UpperPerc; | |
| 692 | 18201396 | qLower = UpperPerc - LowerLoss - LowerEvap - GWFlow; | |
| 693 | |||
| 694 | // --- d(upper zone moisture)/dt = (net upper zone flow) / | ||
| 695 | // (upper zone depth) | ||
| 696 | 18201396 | denom = TotalDepth - x[LOWERDEPTH]; | |
| 697 |
1/2✓ Branch 0 taken 18201396 times.
✗ Branch 1 not taken.
|
18201396 | if (denom > 0.0) |
| 698 | 18201396 | dxdt[THETA] = qUpper / denom; | |
| 699 | else | ||
| 700 | ✗ | dxdt[THETA] = 0.0; | |
| 701 | |||
| 702 | // --- d(lower zone depth)/dt = (net lower zone flow) / | ||
| 703 | // (upper zone moisture deficit) | ||
| 704 | 18201396 | denom = A.porosity - x[THETA]; | |
| 705 |
1/2✓ Branch 0 taken 18201396 times.
✗ Branch 1 not taken.
|
18201396 | if (denom > 0.0) |
| 706 | 18201396 | dxdt[LOWERDEPTH] = qLower / denom; | |
| 707 | else | ||
| 708 | ✗ | dxdt[LOWERDEPTH] = 0.0; | |
| 709 | 18201396 | } | |
| 710 | |||
| 711 | //============================================================================= | ||
| 712 | |||
| 713 | 21234942 | void getEvapRates(double theta, double upperDepth) | |
| 714 | // | ||
| 715 | // Input: theta = moisture content of upper zone | ||
| 716 | // upperDepth = depth of upper zone (ft) | ||
| 717 | // Output: none | ||
| 718 | // Purpose: computes evapotranspiration out of upper & lower zones. | ||
| 719 | // | ||
| 720 | { | ||
| 721 | int p, month; | ||
| 722 | double f; | ||
| 723 | double lowerFrac, upperFrac; | ||
| 724 | |||
| 725 | // --- no GW evaporation when infiltration is occurring | ||
| 726 | 21234942 | UpperEvap = 0.0; | |
| 727 | 21234942 | LowerEvap = 0.0; | |
| 728 |
2/2✓ Branch 0 taken 5650044 times.
✓ Branch 1 taken 15584898 times.
|
21234942 | if ( Infil > 0.0 ) return; |
| 729 | |||
| 730 | // --- get monthly-adjusted upper zone evap fraction | ||
| 731 | 15584898 | upperFrac = A.upperEvapFrac; | |
| 732 | 15584898 | f = 1.0; | |
| 733 | 15584898 | p = A.upperEvapPat; | |
| 734 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15584898 times.
|
15584898 | if ( p >= 0 ) |
| 735 | { | ||
| 736 | ✗ | month = datetime_monthOfYear(getDateTime(NewRunoffTime)); | |
| 737 | ✗ | f = Pattern[p].factor[month-1]; | |
| 738 | } | ||
| 739 | 15584898 | upperFrac *= f; | |
| 740 | |||
| 741 | // --- upper zone evaporation requires that soil moisture | ||
| 742 | // be above the wilting point | ||
| 743 |
2/2✓ Branch 0 taken 15568230 times.
✓ Branch 1 taken 16668 times.
|
15584898 | if ( theta > A.wiltingPoint ) |
| 744 | { | ||
| 745 | // --- actual evap is upper zone fraction applied to max. potential | ||
| 746 | // rate, limited by the available rate after any surface evap | ||
| 747 | 15568230 | UpperEvap = upperFrac * MaxEvap; | |
| 748 |
2/2✓ Branch 0 taken 6855421 times.
✓ Branch 1 taken 8712809 times.
|
15568230 | UpperEvap = MIN(UpperEvap, AvailEvap); |
| 749 | } | ||
| 750 | |||
| 751 | // --- check if lower zone evaporation is possible | ||
| 752 |
1/2✓ Branch 0 taken 15584898 times.
✗ Branch 1 not taken.
|
15584898 | if ( A.lowerEvapDepth > 0.0 ) |
| 753 | { | ||
| 754 | // --- find the fraction of the lower evaporation depth that | ||
| 755 | // extends into the saturated lower zone | ||
| 756 | 15584898 | lowerFrac = (A.lowerEvapDepth - upperDepth) / A.lowerEvapDepth; | |
| 757 |
2/2✓ Branch 0 taken 7517042 times.
✓ Branch 1 taken 8067856 times.
|
15584898 | lowerFrac = MAX(0.0, lowerFrac); |
| 758 |
1/2✓ Branch 0 taken 15584898 times.
✗ Branch 1 not taken.
|
15584898 | lowerFrac = MIN(lowerFrac, 1.0); |
| 759 | |||
| 760 | // --- make the lower zone evap rate proportional to this fraction | ||
| 761 | // and the evap not used in the upper zone | ||
| 762 | 15584898 | LowerEvap = lowerFrac * (1.0 - upperFrac) * MaxEvap; | |
| 763 |
2/2✓ Branch 0 taken 10946936 times.
✓ Branch 1 taken 4637962 times.
|
15584898 | LowerEvap = MIN(LowerEvap, (AvailEvap - UpperEvap)); |
| 764 | } | ||
| 765 | } | ||
| 766 | |||
| 767 | //============================================================================= | ||
| 768 | |||
| 769 | 21234942 | double getUpperPerc(double theta, double upperDepth) | |
| 770 | // | ||
| 771 | // Input: theta = moisture content of upper zone | ||
| 772 | // upperDepth = depth of upper zone (ft) | ||
| 773 | // Output: returns percolation rate (ft/sec) | ||
| 774 | // Purpose: finds percolation rate from upper to lower zone. | ||
| 775 | // | ||
| 776 | { | ||
| 777 | double delta; // unfilled water content of upper zone | ||
| 778 | double dhdz; // avg. change in head with depth | ||
| 779 | double hydcon; // unsaturated hydraulic conductivity | ||
| 780 | |||
| 781 | // --- no perc. from upper zone if no depth or moisture content too low | ||
| 782 |
3/4✓ Branch 0 taken 21234942 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16832980 times.
✓ Branch 3 taken 4401962 times.
|
21234942 | if ( upperDepth <= 0.0 || theta <= A.fieldCapacity ) return 0.0; |
| 783 | |||
| 784 | // --- compute hyd. conductivity as function of moisture content | ||
| 785 | 4401962 | delta = theta - A.porosity; | |
| 786 | 4401962 | hydcon = A.conductivity * exp(delta * A.conductSlope); | |
| 787 | |||
| 788 | // --- compute integral of dh/dz term | ||
| 789 | 4401962 | delta = theta - A.fieldCapacity; | |
| 790 | 4401962 | dhdz = 1.0 + A.tensionSlope * 2.0 * delta / upperDepth; | |
| 791 | |||
| 792 | // --- compute upper zone percolation rate | ||
| 793 | 4401962 | HydCon = hydcon; | |
| 794 | 4401962 | return hydcon * dhdz; | |
| 795 | } | ||
| 796 | |||
| 797 | //============================================================================= | ||
| 798 | |||
| 799 | 21234942 | double getGWFlow(double lowerDepth) | |
| 800 | // | ||
| 801 | // Input: lowerDepth = depth of lower zone (ft) | ||
| 802 | // Output: returns groundwater flow rate (ft/sec) | ||
| 803 | // Purpose: finds groundwater outflow from lower saturated zone. | ||
| 804 | // | ||
| 805 | { | ||
| 806 | double q, t1, t2, t3; | ||
| 807 | |||
| 808 | // --- water table must be above Hstar for flow to occur | ||
| 809 |
2/2✓ Branch 0 taken 267363 times.
✓ Branch 1 taken 20967579 times.
|
21234942 | if ( lowerDepth <= Hstar ) return 0.0; |
| 810 | |||
| 811 | // --- compute groundwater component of flow | ||
| 812 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20967579 times.
|
20967579 | if ( GW->b1 == 0.0 ) t1 = GW->a1; |
| 813 | 20967579 | else t1 = GW->a1 * pow( (lowerDepth - Hstar)*UCF(LENGTH), GW->b1); | |
| 814 | |||
| 815 | // --- compute surface water component of flow | ||
| 816 |
2/2✓ Branch 0 taken 20846505 times.
✓ Branch 1 taken 121074 times.
|
20967579 | if ( GW->b2 == 0.0 ) t2 = GW->a2; |
| 817 |
2/2✓ Branch 0 taken 120166 times.
✓ Branch 1 taken 908 times.
|
121074 | else if (Hsw > Hstar) |
| 818 | { | ||
| 819 | 120166 | t2 = GW->a2 * pow( (Hsw - Hstar)*UCF(LENGTH), GW->b2); | |
| 820 | } | ||
| 821 | 908 | else t2 = 0.0; | |
| 822 | |||
| 823 | // --- compute groundwater/surface water interaction term | ||
| 824 | 20967579 | t3 = GW->a3 * lowerDepth * Hsw * UCF(LENGTH) * UCF(LENGTH); | |
| 825 | |||
| 826 | // --- compute total groundwater flow | ||
| 827 | 20967579 | q = (t1 - t2 + t3) / UCF(GWFLOW); | |
| 828 |
3/4✓ Branch 0 taken 111312 times.
✓ Branch 1 taken 20856267 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 111312 times.
|
20967579 | if ( q < 0.0 && GW->a3 != 0.0 ) q = 0.0; |
| 829 | 20967579 | return q; | |
| 830 | } | ||
| 831 | |||
| 832 | //============================================================================= | ||
| 833 | |||
| 834 | 24 | int getVariableIndex(char* s) | |
| 835 | // | ||
| 836 | // Input: s = name of a groundwater variable | ||
| 837 | // Output: returns index of groundwater variable | ||
| 838 | // Purpose: finds position of GW variable in list of GW variable names. | ||
| 839 | // | ||
| 840 | { | ||
| 841 | int k; | ||
| 842 | |||
| 843 | 24 | k = findmatch(s, GWVarWords); | |
| 844 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | if ( k >= 0 ) return k; |
| 845 | ✗ | return -1; | |
| 846 | } | ||
| 847 | |||
| 848 | //============================================================================= | ||
| 849 | |||
| 850 | 484200 | double getVariableValue(int varIndex) | |
| 851 | // | ||
| 852 | // Input: varIndex = index of a GW variable | ||
| 853 | // Output: returns current value of GW variable | ||
| 854 | // Purpose: finds current value of a GW variable. | ||
| 855 | // | ||
| 856 | { | ||
| 857 |
3/12✓ Branch 0 taken 161400 times.
✓ Branch 1 taken 161400 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 161400 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
484200 | switch (varIndex) |
| 858 | { | ||
| 859 | 161400 | case gwvHGW: return Hgw * UCF(LENGTH); | |
| 860 | 161400 | case gwvHSW: return Hsw * UCF(LENGTH); | |
| 861 | ✗ | case gwvHCB: return Hstar * UCF(LENGTH); | |
| 862 | ✗ | case gwvHGS: return TotalDepth * UCF(LENGTH); | |
| 863 | 161400 | case gwvKS: return A.conductivity * UCF(RAINFALL); | |
| 864 | ✗ | case gwvK: return HydCon * UCF(RAINFALL); | |
| 865 | ✗ | case gwvTHETA:return Theta; | |
| 866 | ✗ | case gwvPHI: return A.porosity; | |
| 867 | ✗ | case gwvFI: return Infil * UCF(RAINFALL); | |
| 868 | ✗ | case gwvFU: return UpperPerc * UCF(RAINFALL); | |
| 869 | ✗ | case gwvA: return Area * UCF(LANDAREA); | |
| 870 | ✗ | default: return 0.0; | |
| 871 | } | ||
| 872 | } | ||
| 873 |