GCC Code Coverage Report

Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 93.4% 454 / 0 / 486
Functions: 100.0% 35 / 0 / 35
Branches: 74.0% 259 / 0 / 350
rdii.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // rdii.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 11/01/21 (Build 5.2.0)
7 // Author: L. Rossman
8 // R. Dickinson (CDM)
9 //
10 // RDII processing functions.
11 //
12 // Note: RDII means rainfall dependent infiltration/inflow,
13 // UH means unit hydrograph.
14 //
15 // Update History
16 // ==============
17 // Build 5.1.007:
18 // - Ignore RDII option implemented.
19 // - Rainfall climate adjustment implemented.
20 // Build 5.1.014:
21 // - Fixes bug related to isUsed property of a unit hydrograph's rain gage.
22 //-----------------------------------------------------------------------------
23 #define _CRT_SECURE_NO_DEPRECATE
24
25 #include <math.h>
26 #include <string.h>
27 #include <stdlib.h>
28 #include "headers.h"
29
30 //-----------------------------------------------------------------------------
31 // Definition of 4-byte integer, 4-byte real and 8-byte real types
32 //-----------------------------------------------------------------------------
33 #define INT4 int
34 #define REAL4 float
35 #define REAL8 double
36 #define FILE_STAMP "SWMM5-RDII"
37
38 //-----------------------------------------------------------------------------
39 // Constants
40 //-----------------------------------------------------------------------------
41 const double ZERO_RDII = 0.0001; // Minimum non-zero RDII inflow (cfs)
42 const char FileStamp[] = FILE_STAMP;
43
44 //-----------------------------------------------------------------------------
45 // Data Structures
46 //-----------------------------------------------------------------------------
47 enum FileTypes {BINARY, TEXT}; // File mode types
48
49 typedef struct // Data for a single unit hydrograph
50 { // -------------------------------------
51 double* pastRain; // array of past rainfall values
52 char* pastMonth; // month in which past rainfall occurred
53 int period; // current UH time period
54 int hasPastRain; // true if > 0 past periods with rain
55 int maxPeriods; // max. past rainfall periods
56 long drySeconds; // time since last nonzero rainfall
57 double iaUsed; // initial abstraction used (in or mm)
58 } TUHData;
59
60 typedef struct // Data for a unit hydrograph group
61 { //---------------------------------
62 int isUsed; // true if UH group used by any nodes
63 int rainInterval; // time interval for RDII processing (sec)
64 double area; // sewered area covered by UH's gage (ft2)
65 double rdii; // rdii flow (in rainfall units)
66 DateTime gageDate; // calendar date of rain gage period
67 DateTime lastDate; // date of last rdii computed
68 TUHData uh[3]; // data for each unit hydrograph
69 } TUHGroup;
70
71 //-----------------------------------------------------------------------------
72 // Shared Variables
73 //-----------------------------------------------------------------------------
74 static TUHGroup* UHGroup; // processing data for each UH group
75 static int RdiiStep; // RDII time step (sec)
76 static int NumRdiiNodes; // number of nodes w/ RDII data
77 static int* RdiiNodeIndex; // indexes of nodes w/ RDII data
78 static REAL4* RdiiNodeFlow; // inflows for nodes with RDII
79 static int RdiiFlowUnits; // RDII flow units code
80 static DateTime RdiiStartDate; // start date of RDII inflow period
81 static DateTime RdiiEndDate; // end date of RDII inflow period
82 static double TotalRainVol; // total rainfall volume (ft3)
83 static double TotalRdiiVol; // total RDII volume (ft3)
84 static int RdiiFileType; // type (binary/text) of RDII file
85
86 //-----------------------------------------------------------------------------
87 // Imported Variables
88 //-----------------------------------------------------------------------------
89 #ifdef __cplusplus
90 extern const double Qcf[]; // flow units conversion factors
91 // (see swmm5.c)
92 #else
93 extern double Qcf[]; // flow units conversion factors
94 // (see swmm5.c)
95 #endif
96
97 //-----------------------------------------------------------------------------
98 // External functions (declared in funcs.h)
99 //-----------------------------------------------------------------------------
100 // rdii_readRdiiInflow (called from parseLine in input.c)
101 // rdii_deleteRdiiInflow (called from deleteObjects in project.c)
102 // rdii_initUnitHyd (called from createObjects in project.c)
103 // rdii_readUnitHydParams (called from parseLine in input.c)
104 // rdii_openRdii (called from rain_open)
105 // rdii_closeRdii (called from rain_close)
106 // rdii_getNumRdiiFlows (called from addRdiiInflows in routing.c)
107 // rdii_getRdiiFlow (called from addRdiiInflows in routing.c)
108
109 //-----------------------------------------------------------------------------
110 // Function Declarations
111 //-----------------------------------------------------------------------------
112 // --- functions used to create a RDII file
113 static int readOldUHFormat(int j, int m, char* tok[], int ntoks);
114 static void setUnitHydParams(int j, int i, int m, double x[]);
115 static void createRdiiFile(void);
116 static int getNumRdiiNodes(void);
117 static void validateRdii(void);
118
119 static void openRdiiProcessor(void);
120 static int allocRdiiMemory(void);
121 static int getRainInterval(int i);
122 static int getMaxPeriods(int i, int k);
123 static void initGageData(void);
124 static void initUnitHydData(void);
125 static int openNewRdiiFile(void);
126 static void getRainfall(DateTime currentDate);
127
128 static double applyIA(int j, int k, DateTime aDate, double dt,
129 double rainDepth);
130 static void updateDryPeriod(int j, int k, double rain, int gageInterval);
131 static void getUnitHydRdii(DateTime currentDate);
132 static double getUnitHydConvol(int j, int k, int gageInterval);
133 static double getUnitHydOrd(int j, int m, int k, double t);
134
135 static int getNodeRdii(void);
136 static void saveRdiiFlows(DateTime currentDate);
137 static void closeRdiiProcessor(void);
138 static void freeRdiiMemory(void);
139
140 // --- functions used to read an existing RDII file
141 static int readRdiiFileHeader(void);
142 static void readRdiiFlows(void);
143
144 static void openRdiiTextFile(void);
145 static int readRdiiTextFileHeader(void);
146 static void readRdiiTextFlows(void);
147
148 //=============================================================================
149 // Management of RDII-Related Data
150 //=============================================================================
151
152 3046 int rdii_readRdiiInflow(char* tok[], int ntoks)
153 //
154 // Input: tok[] = array of string tokens
155 // ntoks = number of tokens
156 // Output: returns an error code
157 // Purpose: reads properties of an RDII inflow from a line of input.
158 //
159 {
160 int j, k;
161 double a;
162 TRdiiInflow* inflow;
163
164 // --- check for proper number of items
165
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3046 times.
 3046 if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
166
167 // --- check that node receiving RDII exists
168 3046 j = project_findObject(NODE, tok[0]);
169
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3046 times.
 3046 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
170
171 // --- check that RDII unit hydrograph exists
172 3046 k = project_findObject(UNITHYD, tok[1]);
173
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3046 times.
 3046 if ( k < 0 ) return error_setInpError(ERR_NAME, tok[1]);
174
175 // --- read in sewer area value
176
2/4
✓ Branch 1 taken 3046 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 3046 times.
 3046 if ( !getDouble(tok[2], &a) || a < 0.0 )
177 ✗ return error_setInpError(ERR_NUMBER, tok[2]);
178
179 // --- create the RDII inflow object if it doesn't already exist
180 3046 inflow = Node[j].rdiiInflow;
181
1/2
✓ Branch 0 taken 3046 times.
✗ Branch 1 not taken.
 3046 if ( inflow == NULL )
182 {
183 3046 inflow = (TRdiiInflow *) malloc(sizeof(TRdiiInflow));
184
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3046 times.
 3046 if ( !inflow ) return error_setInpError(ERR_MEMORY, "");
185 }
186
187 // --- assign UH & area to inflow object
188 3046 inflow->unitHyd = k;
189 3046 inflow->area = a / UCF(LANDAREA);
190
191 // --- assign inflow object to node
192 3046 Node[j].rdiiInflow = inflow;
193 3046 return 0;
194 }
195
196 //=============================================================================
197
198 220 void rdii_initUnitHyd(int j)
199 //
200 // Input: j = UH group index
201 // Output: none
202 // Purpose: initializes properties of a unit hydrograph group.
203 //
204 {
205 int i; // individual UH index
206 int m; // month index
207
208
2/2
✓ Branch 0 taken 2640 times.
✓ Branch 1 taken 220 times.
 2860 for ( m=0; m<12; m++)
209 {
210
2/2
✓ Branch 0 taken 7920 times.
✓ Branch 1 taken 2640 times.
 10560 for (i=0; i<3; i++)
211 {
212 7920 UnitHyd[j].iaMax[m][i] = 0.0;
213 7920 UnitHyd[j].iaRecov[m][i] = 0.0;
214 7920 UnitHyd[j].iaInit[m][i] = 0.0;
215 7920 UnitHyd[j].r[m][i] = 0.0;
216 7920 UnitHyd[j].tPeak[m][i] = 0;
217 7920 UnitHyd[j].tBase[m][i] = 0;
218 }
219 }
220 220 }
221
222 //=============================================================================
223
224 8006 int rdii_readUnitHydParams(char* tok[], int ntoks)
225 //
226 // Input: tok[] = array of string tokens
227 // ntoks = number of tokens
228 // Output: returns an error code
229 // Purpose: reads parameters of an RDII unit hydrograph from a line of input.
230 //
231 {
232 int i, j, k, m, g;
233 double x[6];
234
235 // --- check that RDII UH object exists in database
236 8006 j = project_findObject(UNITHYD, tok[0]);
237
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8006 times.
 8006 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
238
239 // --- assign UH ID to name in hash table
240
2/2
✓ Branch 0 taken 220 times.
✓ Branch 1 taken 7786 times.
 8006 if ( UnitHyd[j].ID == NULL )
241 220 UnitHyd[j].ID = project_findID(UNITHYD, tok[0]);
242
243 // --- line has 2 tokens; assign rain gage to UH object
244
2/2
✓ Branch 0 taken 220 times.
✓ Branch 1 taken 7786 times.
 8006 if ( ntoks == 2 )
245 {
246 220 g = project_findObject(GAGE, tok[1]);
247
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 220 times.
 220 if ( g < 0 ) return error_setInpError(ERR_NAME, tok[1]);
248 220 UnitHyd[j].rainGage = g;
249 220 Gage[g].isUsed = TRUE;
250 220 return 0;
251 }
252
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7786 times.
 7786 else if ( ntoks < 6 ) return error_setInpError(ERR_ITEMS, "");
253
254 // --- find which month UH params apply to
255 7786 m = datetime_findMonth(tok[1]);
256
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 7774 times.
 7786 if ( m == 0 )
257 {
258
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 if ( !match(tok[1], w_ALL) )
259 ✗ return error_setInpError(ERR_KEYWORD, tok[1]);
260 }
261
262 // --- find type of UH being specified
263 7786 k = findmatch(tok[2], UHTypeWords);
264
265 // --- if no type match, try using older UH line format
266
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 7785 times.
 7786 if ( k < 0 ) return readOldUHFormat(j, m, tok, ntoks);
267
268 // --- read the R-T-K parameters
269
2/2
✓ Branch 0 taken 23355 times.
✓ Branch 1 taken 7785 times.
 31140 for ( i = 0; i < 3; i++ )
270 {
271
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 23355 times.
 23355 if ( ! getDouble(tok[i+3], &x[i]) )
272 ✗ return error_setInpError(ERR_NUMBER, tok[i+3]);
273 }
274
275 // --- read the IA parameters if present
276
2/2
✓ Branch 0 taken 23355 times.
✓ Branch 1 taken 7785 times.
 31140 for (i = 3; i < 6; i++)
277 {
278 23355 x[i] = 0.0;
279
1/2
✓ Branch 0 taken 23355 times.
✗ Branch 1 not taken.
 23355 if ( ntoks > i+3 )
280 {
281
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 23355 times.
 23355 if ( ! getDouble(tok[i+3], &x[i]) )
282 ✗ return error_setInpError(ERR_NUMBER, tok[i+2]);
283 }
284 }
285
286 // --- save UH params
287 7785 setUnitHydParams(j, k, m, x);
288 7785 return 0;
289 }
290
291 //=============================================================================
292
293 1 int readOldUHFormat(int j, int m, char* tok[], int ntoks)
294 //
295 // Input: j = unit hydrograph index
296 // m = month of year (0 = all months)
297 // tok[] = array of string tokens
298 // ntoks = number of tokens
299 // Output: returns an error code
300 // Purpose: reads parameters of a set of RDII unit hydrographs from a line of
301 // input.
302 //
303 {
304 int i, k;
305 double p[9], x[6];
306
307 // --- check for proper number of tokens
308
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( ntoks < 11 ) return error_setInpError(ERR_ITEMS, "");
309
310 // --- read 3 sets of r-t-k values
311
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 1 time.
 10 for ( i = 0; i < 9; i++ )
312 {
313
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
 9 if ( ! getDouble(tok[i+2], &p[i]) )
314 ✗ return error_setInpError(ERR_NUMBER, tok[i+2]);
315 }
316
317 // --- read initial abstraction parameters
318
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
 4 for (i = 0; i < 3; i++)
319 {
320 3 x[i+3] = 0.0;
321
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if ( ntoks > i+11 )
322 {
323
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
 3 if ( ! getDouble(tok[i+11], &x[i+3]) )
324 ✗ return error_setInpError(ERR_NUMBER, tok[i+11]);
325 }
326 }
327
328 // --- save UH parameters
329
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
 4 for ( k = 0; k < 3; k++)
330 {
331
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
 12 for ( i = 0; i < 3; i++)
332 {
333 9 x[i] = p[3*k + i];
334 9 setUnitHydParams(j, k, m, x);
335 }
336 }
337 1 return 0;
338 }
339
340 //=============================================================================
341
342 7794 void setUnitHydParams(int j, int i, int m, double x[])
343 //
344 // Input: j = unit hydrograph index
345 // i = type of UH response (short, medium or long term)
346 // m = month of year (0 = all months)
347 // x = array of UH parameters
348 // Output: none
349 // Purpose: assigns parameters to a unit hydrograph for a specified month of year.
350 //
351 {
352 int m1, m2; // start/end month indexes
353 double t, // UH time to peak (hrs)
354 k, // UH k-value
355 tBase; // UH base time (hrs)
356
357 // --- find range of months that share same parameter values
358
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 7782 times.
 7794 if ( m == 0 )
359 {
360 12 m1 = 0;
361 12 m2 = 11;
362 }
363 else
364 {
365 7782 m1 = m-1;
366 7782 m2 = m1;
367 }
368
369 // --- for each month in the range
370
2/2
✓ Branch 0 taken 7926 times.
✓ Branch 1 taken 7794 times.
 15720 for (m=m1; m<=m2; m++)
371 {
372 // --- set UH response ratio, time to peak, & base time
373 7926 UnitHyd[j].r[m][i] = x[0];
374 7926 t = x[1];
375 7926 k = x[2];
376 7926 tBase = t * (1.0 + k); // hours
377 7926 UnitHyd[j].tPeak[m][i] = (long)(t * 3600.); // seconds
378 7926 UnitHyd[j].tBase[m][i] = (long)(tBase * 3600.); // seconds
379
380 // -- set initial abstraction parameters
381 7926 UnitHyd[j].iaMax[m][i] = x[3];
382 7926 UnitHyd[j].iaRecov[m][i] = x[4];
383 7926 UnitHyd[j].iaInit[m][i] = x[5];
384 }
385 7794 }
386
387 //=============================================================================
388
389 10183 void rdii_deleteRdiiInflow(int j)
390 //
391 // Input: j = node index
392 // Output: none
393 // Purpose: deletes the RDII inflow object for a node.
394 //
395 {
396
2/2
✓ Branch 0 taken 3046 times.
✓ Branch 1 taken 7137 times.
 10183 if ( Node[j].rdiiInflow )
397 {
398 3046 free(Node[j].rdiiInflow);
399 3046 Node[j].rdiiInflow = NULL;
400 }
401 10183 }
402
403
404 //=============================================================================
405 // Reading Inflow Data From a RDII File
406 //=============================================================================
407
408 58 void rdii_openRdii()
409 //
410 // Input: none
411 // Output: none
412 // Purpose: opens an exisiting RDII interface file or creates a new one.
413 //
414 {
415 58 char fStamp[] = FILE_STAMP;
416
417 58 RdiiNodeIndex = NULL;
418 58 RdiiNodeFlow = NULL;
419 58 NumRdiiNodes = 0;
420 58 RdiiStartDate = NO_DATE;
421
422 // --- create the RDII file if existing file not being used
423
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 58 times.
 111 if ( IgnoreRDII ) return;
424
2/2
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 2 times.
 58 if ( Frdii.mode != USE_FILE ) createRdiiFile();
425
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
 58 if ( Frdii.mode == NO_FILE || ErrorCode ) return;
426
427 // --- try to open the RDII file in binary mode
428 5 Frdii.file = fopen(Frdii.name, "rb");
429
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( Frdii.file == NULL)
430 {
431 ✗ if ( Frdii.mode == SCRATCH_FILE )
432 {
433 ✗ report_writeErrorMsg(ERR_RDII_FILE_SCRATCH, "");
434 }
435 else
436 {
437 ✗ report_writeErrorMsg(ERR_RDII_FILE_OPEN, Frdii.name);
438 }
439 ✗ return;
440 }
441
442 // --- check for valid file stamp
443 5 fread(fStamp, sizeof(char), strlen(FileStamp), Frdii.file);
444
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 time.
 5 if ( strcmp(fStamp, FileStamp) == 0 )
445 {
446 4 RdiiFileType = BINARY;
447 4 ErrorCode = readRdiiFileHeader();
448 }
449
450 // --- if stamp invalid try to open the file in text mode
451 else
452 {
453 1 fclose(Frdii.file);
454 1 RdiiFileType = TEXT;
455 1 openRdiiTextFile();
456 }
457
458 // --- catch any error
459
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( ErrorCode )
460 {
461 ✗ report_writeErrorMsg(ErrorCode, Frdii.name);
462 }
463
464 // --- read the first set of RDII flows form the file
465 5 else readRdiiFlows();
466 }
467
468 //=============================================================================
469
470 1 void openRdiiTextFile()
471 {
472 // --- try to open the RDII file in text mode
473 1 Frdii.file = fopen(Frdii.name, "rt");
474
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( Frdii.file == NULL)
475 {
476 ✗ if ( Frdii.mode == SCRATCH_FILE )
477 {
478 ✗ report_writeErrorMsg(ERR_RDII_FILE_SCRATCH, "");
479 }
480 else
481 {
482 ✗ report_writeErrorMsg(ERR_RDII_FILE_OPEN, Frdii.name);
483 }
484 ✗ return;
485 }
486
487 // --- read header records from file
488 1 ErrorCode = readRdiiTextFileHeader();
489
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( ErrorCode )
490 {
491 ✗ report_writeErrorMsg(ErrorCode, Frdii.name);
492 }
493 }
494
495 //=============================================================================
496
497 58 void rdii_closeRdii()
498 //
499 // Input: none
500 // Output: none
501 // Purpose: closes the RDII interface file.
502 //
503 {
504
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 53 times.
 58 if ( Frdii.file ) fclose(Frdii.file);
505
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 56 times.
 58 if ( Frdii.mode == SCRATCH_FILE ) remove(Frdii.name);
506
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 53 times.
 58 FREE(RdiiNodeIndex);
507
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 53 times.
 58 FREE(RdiiNodeFlow);
508 58 }
509
510 //=============================================================================
511
512 1030493 int rdii_getNumRdiiFlows(DateTime aDate)
513 //
514 // Input: aDate = current date/time
515 // Output: returns 0 if no RDII flow or number of nodes with RDII inflows
516 // Purpose: finds number of RDII inflows at a specified date.
517 //
518 {
519 // --- default result is 0 indicating no RDII inflow at specified date
520
2/2
✓ Branch 0 taken 1006205 times.
✓ Branch 1 taken 24288 times.
 1030493 if ( NumRdiiNodes == 0 ) return 0;
521
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24288 times.
 24288 if ( !Frdii.file ) return 0;
522
523 // --- keep reading RDII file as need be
524
2/2
✓ Branch 1 taken 22138 times.
✓ Branch 2 taken 3308 times.
 25446 while ( !feof(Frdii.file) )
525 {
526 // --- return if date of current RDII inflow not reached yet
527
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 22138 times.
 22138 if ( RdiiStartDate == NO_DATE ) return 0;
528
2/2
✓ Branch 0 taken 7418 times.
✓ Branch 1 taken 14720 times.
 22138 if ( aDate < RdiiStartDate ) return 0;
529
530 // --- return RDII node count if specified date falls
531 // within time interval of current RDII inflow
532
2/2
✓ Branch 0 taken 13562 times.
✓ Branch 1 taken 1158 times.
 14720 if ( aDate < RdiiEndDate ) return NumRdiiNodes;
533
534 // --- otherwise get next date and RDII flow values from file
535 1158 else readRdiiFlows();
536 }
537 3308 return 0;
538 }
539
540 //=============================================================================
541
542 10971860 void rdii_getRdiiFlow(int i, int* j, double* q)
543 //
544 // Input: i = RDII node index
545 // j = pointer to project node index
546 // q = pointer to RDII flow rate
547 // Output: sets node index and RDII inflow for node
548 // Purpose: finds index and current RDII inflow for an RDII node.
549 //
550 {
551
2/4
✓ Branch 0 taken 10971860 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10971860 times.
✗ Branch 3 not taken.
 10971860 if ( i >= 0 && i < NumRdiiNodes )
552 {
553 10971860 *j = RdiiNodeIndex[i];
554 10971860 *q = RdiiNodeFlow[i];
555 }
556 10971860 }
557
558 //=============================================================================
559
560 4 int readRdiiFileHeader()
561 //
562 // Input: none
563 // Output: returns error code
564 // Purpose: reads header information from a binary RDII file.
565 //
566 {
567 int i, j;
568
569 // --- extract time step and number of RDII nodes
570 4 fread(&RdiiStep, sizeof(INT4), 1, Frdii.file);
571
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( RdiiStep <= 0 ) return ERR_RDII_FILE_FORMAT;
572 4 fread(&NumRdiiNodes, sizeof(INT4), 1, Frdii.file);
573
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( NumRdiiNodes <= 0 ) return ERR_RDII_FILE_FORMAT;
574
575 // --- allocate memory for RdiiNodeIndex & RdiiNodeFlow arrays
576 4 RdiiNodeIndex = (int *) calloc(NumRdiiNodes, sizeof(int));
577
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( !RdiiNodeIndex ) return ERR_MEMORY;
578 4 RdiiNodeFlow = (REAL4 *) calloc(NumRdiiNodes, sizeof(REAL4));
579
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( !RdiiNodeFlow ) return ERR_MEMORY;
580
581 // --- read indexes of RDII nodes
582
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( feof(Frdii.file) ) return ERR_RDII_FILE_FORMAT;
583 4 fread(RdiiNodeIndex, sizeof(INT4), NumRdiiNodes, Frdii.file);
584
2/2
✓ Branch 0 taken 3046 times.
✓ Branch 1 taken 4 times.
 3050 for ( i=0; i<NumRdiiNodes; i++ )
585 {
586 3046 j = RdiiNodeIndex[i];
587
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3046 times.
 3046 if ( Node[j].rdiiInflow == NULL ) return ERR_RDII_FILE_FORMAT;
588 }
589
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( feof(Frdii.file) ) return ERR_RDII_FILE_FORMAT;
590 4 return 0;
591 }
592
593 //=============================================================================
594
595 1 int readRdiiTextFileHeader()
596 //
597 // Input: none
598 // Output: returns error code
599 // Purpose: reads header information from a text RDII file.
600 //
601 {
602 int i, j;
603 char line[MAXLINE+1]; // line from RDII data file
604 1 char s1[MAXLINE+1] = ""; // general string variable
605 char s2[MAXLINE+1];
606
607 // --- check for correct file type
608 1 fgets(line, MAXLINE, Frdii.file);
609
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
 1 if ( !sscanf(line, "%s", s1) || strcmp(s1, "SWMM5") != 0 )
610 ✗ return ERR_RDII_FILE_FORMAT;
611
612 // --- skip title line
613 1 fgets(line, MAXLINE, Frdii.file);
614
615 // --- read RDII UH time step interval (sec)
616 1 RdiiStep = 0;
617 1 fgets(line, MAXLINE, Frdii.file);
618
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
 1 if ( !sscanf(line, "%d", &RdiiStep) || RdiiStep <= 0 )
619 ✗ return ERR_RDII_FILE_FORMAT;
620
621 // --- skip over line with number of constituents (= 1 for RDII)
622 1 fgets(line, MAXLINE, Frdii.file);
623
624 // --- read flow units
625 1 fgets(line, MAXLINE, Frdii.file);
626
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( sscanf(line, "%s %s", s1, s2) < 2 )
627 ✗ return ERR_RDII_FILE_FORMAT;
628 1 RdiiFlowUnits = findmatch(s2, FlowUnitWords);
629
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( RdiiFlowUnits < 0 ) return ERR_RDII_FILE_FORMAT;
630
631 // --- read number of RDII nodes
632 1 fgets(line, MAXLINE, Frdii.file);
633
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
 1 if ( sscanf(line, "%d", &NumRdiiNodes) < 1 || NumRdiiNodes <= 0 )
634 ✗ return ERR_RDII_FILE_FORMAT;
635
636 // --- allocate memory for RdiiNodeIndex & RdiiNodeFlow arrays
637 1 RdiiNodeIndex = (int *) calloc(NumRdiiNodes, sizeof(int));
638
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( !RdiiNodeIndex ) return ERR_MEMORY;
639 1 RdiiNodeFlow = (REAL4 *) calloc(NumRdiiNodes, sizeof(REAL4));
640
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( !RdiiNodeFlow ) return ERR_MEMORY;
641
642 // --- read names of RDII nodes from file & save their indexes
643
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
 2 for ( i=0; i<NumRdiiNodes; i++ )
644 {
645
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if ( feof(Frdii.file) ) return ERR_RDII_FILE_FORMAT;
646 1 fgets(line, MAXLINE, Frdii.file);
647
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( !sscanf(line, "%s", s1) )
648 ✗ return ERR_RDII_FILE_FORMAT;
649 1 j = project_findObject(NODE, s1);
650
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( j < 0 )
651 ✗ return ERR_RDII_FILE_FORMAT;
652 1 RdiiNodeIndex[i] = j;
653 }
654
655 // --- skip column heading line
656
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if ( feof(Frdii.file) ) return ERR_RDII_FILE_FORMAT;
657 1 fgets(line, MAXLINE, Frdii.file);
658 1 return 0;
659 }
660
661 //=============================================================================
662
663 1163 void readRdiiFlows()
664 //
665 // Input: none
666 // Output: none
667 // Purpose: reads date and flow values of next RDII inflows from RDII file.
668 //
669 {
670
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1157 times.
 1163 if ( RdiiFileType == TEXT ) readRdiiTextFlows();
671 else
672 {
673 1157 RdiiStartDate = NO_DATE;
674 1157 RdiiEndDate = NO_DATE;
675
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1157 times.
 1157 if ( feof(Frdii.file) ) return;
676 1157 fread(&RdiiStartDate, sizeof(DateTime), 1, Frdii.file);
677
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1156 times.
 1157 if ( RdiiStartDate == NO_DATE ) return;
678 1156 if ( fread(RdiiNodeFlow, sizeof(REAL4), NumRdiiNodes, Frdii.file)
679
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1156 times.
 1156 < (size_t)NumRdiiNodes ) RdiiStartDate = NO_DATE;
680 1156 else RdiiEndDate = datetime_addSeconds(RdiiStartDate, RdiiStep);
681 }
682 }
683
684 //=============================================================================
685
686 6 void readRdiiTextFlows()
687 //
688 // Input: none
689 // Output: none
690 // Purpose: reads date and flow values of next RDII inflows from RDII file.
691 //
692 {
693 int i, n;
694 6 int yr = 0, mon = 0, day = 0,
695 6 hr = 0, min = 0, sec = 0; // year, month, day, hour, minute, second
696 double x; // RDII flow in original units
697 char line[MAXLINE+1]; // line from RDII data file
698 char s[MAXLINE+1]; // node ID label (not used)
699
700 6 RdiiStartDate = NO_DATE;
701
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 6 times.
 12 for (i=0; i<NumRdiiNodes; i++)
702 {
703
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
 6 if ( feof(Frdii.file) ) return;
704 6 fgets(line, MAXLINE, Frdii.file);
705 6 n = sscanf(line, "%s %d %d %d %d %d %d %lf",
706 s, &yr, &mon, &day, &hr, &min, &sec, &x);
707
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 if ( n < 8 ) return;
708 6 RdiiNodeFlow[i] = (REAL4)(x / Qcf[RdiiFlowUnits]);
709 }
710 6 RdiiStartDate = datetime_encodeDate(yr, mon, day) +
711 6 datetime_encodeTime(hr, min, sec);
712 6 RdiiEndDate = datetime_addSeconds(RdiiStartDate, RdiiStep);
713 }
714
715
716 //=============================================================================
717 // Creation of a RDII Interface File
718 //=============================================================================
719
720 56 void createRdiiFile()
721 //
722 // Input: none
723 // Output: none
724 // Purpose: computes time history of RDII inflows and saves them to file.
725 //
726 {
727 int hasRdii; // true when total RDII > 0
728 double elapsedTime; // current elapsed time (sec)
729 double duration; // duration being analyzed (sec)
730 DateTime currentDate; // current calendar date/time
731
732 // --- set RDII reporting time step to Runoff wet step
733 56 RdiiStep = WetStep;
734
735 // --- count nodes with RDII data
736 56 NumRdiiNodes = getNumRdiiNodes();
737
738 // --- if no RDII nodes then re-set RDII file usage to NO_FILE
739
2/2
✓ Branch 0 taken 53 times.
✓ Branch 1 taken 3 times.
 56 if ( NumRdiiNodes == 0 )
740 {
741 53 Frdii.mode = NO_FILE;
742 53 return;
743 }
744
745 // --- otherwise set file usage to SCRATCH if originally set to NO_FILE
746
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
 3 else if ( Frdii.mode == NO_FILE ) Frdii.mode = SCRATCH_FILE;
747
748 // --- validate RDII data
749 3 validateRdii();
750 3 initGageData();
751
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if ( ErrorCode ) return;
752
753 // --- open RDII processing system
754 3 openRdiiProcessor();
755
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if ( !ErrorCode )
756 {
757 // --- initialize rain gage & UH processing data
758 3 initUnitHydData();
759
760 // --- convert total simulation duration from millisec to sec
761 3 duration = TotalDuration / 1000.0;
762
763 // --- examine rainfall record over each RdiiStep time step
764 3 elapsedTime = 0.0;
765
3/4
✓ Branch 0 taken 3039 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3039 times.
✗ Branch 3 not taken.
 3042 while ( elapsedTime <= duration && !ErrorCode )
766 {
767 // --- compute current calendar date/time
768 3039 currentDate = StartDateTime + elapsedTime / SECperDAY;
769
770 // --- update rainfall at all rain gages
771 3039 getRainfall(currentDate);
772
773 // --- compute convolutions of past rainfall with UH's
774 3039 getUnitHydRdii(currentDate);
775
776 // --- find RDII at all nodes
777 3039 hasRdii = getNodeRdii();
778
779 // --- save RDII at all nodes to file for current date
780
2/2
✓ Branch 0 taken 1130 times.
✓ Branch 1 taken 1909 times.
 3039 if ( hasRdii ) saveRdiiFlows(currentDate);
781
782 // --- advance one time step
783 3039 elapsedTime += RdiiStep;
784 }
785 }
786
787 // --- close RDII processing system
788 3 closeRdiiProcessor();
789 }
790
791 //=============================================================================
792
793 56 int getNumRdiiNodes()
794 //
795 // Input: none
796 // Output: returns node count
797 // Purpose: counts number of nodes that receive RDII inflow.
798 //
799 {
800 int j, // node index
801 n; // node count
802
803 56 n = 0;
804
2/2
✓ Branch 0 taken 9030 times.
✓ Branch 1 taken 56 times.
 9086 for (j=0; j<Nobjects[NODE]; j++)
805 {
806
2/2
✓ Branch 0 taken 2031 times.
✓ Branch 1 taken 6999 times.
 9030 if ( Node[j].rdiiInflow ) n++;
807 }
808 56 return n;
809 }
810
811 //=============================================================================
812
813 3 void validateRdii()
814 //
815 // Input: none
816 // Output: none
817 // Purpose: validates UH and RDII inflow object data.
818 //
819 {
820 int i, // node index
821 j, // UH group index
822 k, // individual UH index
823 m; // month index
824 double rsum; // sum of UH r-values
825 // long gageInterval; // rain gage time interval
826
827 // --- check each unit hydrograph for consistency
828
2/2
✓ Branch 0 taken 146 times.
✓ Branch 1 taken 3 times.
 149 for (j=0; j<Nobjects[UNITHYD]; j++)
829 {
830
2/2
✓ Branch 0 taken 1752 times.
✓ Branch 1 taken 146 times.
 1898 for (m=0; m<12; m++)
831 {
832 1752 rsum = 0.0;
833
2/2
✓ Branch 0 taken 5256 times.
✓ Branch 1 taken 1752 times.
 7008 for (k=0; k<3; k++)
834 {
835 // --- if no base time then UH doesn't exist
836
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5256 times.
 5256 if ( UnitHyd[j].tBase[m][k] == 0 ) continue;
837
838 // --- restriction on time to peak being less than the
839 // rain gage's recording interval no longer applies
840
841 // --- can't have negative UH parameters
842
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5256 times.
 5256 if ( UnitHyd[j].tPeak[m][k] < 0.0 )
843 {
844 ✗ report_writeErrorMsg(ERR_UNITHYD_TIMES, UnitHyd[j].ID);
845 }
846
847 // --- can't have negative UH response ratio
848
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5256 times.
 5256 if ( UnitHyd[j].r[m][k] < 0.0 )
849 {
850 ✗ report_writeErrorMsg(ERR_UNITHYD_RATIOS, UnitHyd[j].ID);
851 }
852 5256 else rsum += UnitHyd[j].r[m][k];
853 }
854
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1752 times.
 1752 if ( rsum > 1.01 )
855 {
856 ✗ report_writeErrorMsg(ERR_UNITHYD_RATIOS, UnitHyd[j].ID);
857 }
858 }
859 }
860
861 // --- check each node's RDII inflow object
862
2/2
✓ Branch 0 taken 2308 times.
✓ Branch 1 taken 3 times.
 2311 for (i=0; i<Nobjects[NODE]; i++)
863 {
864
2/2
✓ Branch 0 taken 2031 times.
✓ Branch 1 taken 277 times.
 2308 if ( Node[i].rdiiInflow )
865 {
866 // --- check that sewer area is non-negative
867
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2031 times.
 2031 if ( Node[i].rdiiInflow->area < 0.0 )
868 {
869 ✗ report_writeErrorMsg(ERR_RDII_AREA, Node[i].ID);
870 }
871 }
872 }
873 3 }
874
875 //=============================================================================
876
877 3 void openRdiiProcessor()
878 //
879 // Input: none
880 // Output: none
881 // Purpose: opens RDII processing system.
882 //
883 {
884 int j; // object index
885 int n; // RDII node count
886
887 // --- set RDII processing arrays to NULL
888 3 UHGroup = NULL;
889 3 RdiiNodeIndex = NULL;
890 3 RdiiNodeFlow = NULL;
891 3 TotalRainVol = 0.0;
892 3 TotalRdiiVol = 0.0;
893
894 // --- allocate memory used for RDII processing
895
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
 3 if ( !allocRdiiMemory() )
896 {
897 ✗ report_writeErrorMsg(ERR_MEMORY, "");
898 ✗ return;
899 }
900
901 // --- open & initialize RDII file
902
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
 3 if ( !openNewRdiiFile() )
903 {
904 ✗ report_writeErrorMsg(ERR_RDII_FILE_SCRATCH, "");
905 ✗ return;
906 }
907
908 // --- identify index of each node with RDII inflow
909 3 n = 0;
910
2/2
✓ Branch 0 taken 2308 times.
✓ Branch 1 taken 3 times.
 2311 for (j=0; j<Nobjects[NODE]; j++)
911 {
912
3/4
✓ Branch 0 taken 2031 times.
✓ Branch 1 taken 277 times.
✓ Branch 2 taken 2031 times.
✗ Branch 3 not taken.
 2308 if ( Node[j].rdiiInflow && RdiiNodeIndex != NULL )
913 {
914 2031 RdiiNodeIndex[n] = j;
915 2031 n++;
916 }
917 }
918 }
919
920 //=============================================================================
921
922 3 int allocRdiiMemory()
923 //
924 // Input: none
925 // Output: returns TRUE if successful, FALSE if not
926 // Purpose: allocates memory used for RDII processing .
927 //
928 //
929 {
930 int i; // UH group index
931 int k; // UH index
932 int n; // number of past rain periods
933
934 // --- allocate memory for RDII processing data for UH groups
935 3 UHGroup = (TUHGroup *) calloc(Nobjects[UNITHYD], sizeof(TUHGroup));
936
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if ( !UHGroup ) return FALSE;
937
938 // --- allocate memory for past rainfall data for each UH in each group
939
2/2
✓ Branch 0 taken 146 times.
✓ Branch 1 taken 3 times.
 149 for (i=0; i<Nobjects[UNITHYD]; i++)
940 {
941 146 UHGroup[i].rainInterval = getRainInterval(i);
942
2/2
✓ Branch 0 taken 438 times.
✓ Branch 1 taken 146 times.
 584 for (k=0; k<3; k++)
943 {
944 438 UHGroup[i].uh[k].pastRain = NULL;
945 438 UHGroup[i].uh[k].pastMonth = NULL;
946 438 UHGroup[i].uh[k].maxPeriods = getMaxPeriods(i, k);
947 438 n = UHGroup[i].uh[k].maxPeriods;
948
1/2
✓ Branch 0 taken 438 times.
✗ Branch 1 not taken.
 438 if ( n > 0 )
949 {
950 438 UHGroup[i].uh[k].pastRain =
951 438 (double *) calloc(n, sizeof(double));
952
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 438 times.
 438 if ( !UHGroup[i].uh[k].pastRain ) return FALSE;
953 438 UHGroup[i].uh[k].pastMonth =
954 438 (char *) calloc(n, sizeof(char));
955
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 438 times.
 438 if ( !UHGroup[i].uh[k].pastMonth ) return FALSE;
956 }
957 }
958 }
959
960 // --- allocate memory for RDII indexes & inflow at each node w/ RDII data
961 3 RdiiNodeIndex = (int *) calloc(NumRdiiNodes, sizeof(int));
962
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if ( !RdiiNodeIndex ) return FALSE;
963 3 RdiiNodeFlow = (REAL4 *) calloc(NumRdiiNodes, sizeof(REAL4));
964
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if ( !RdiiNodeFlow ) return FALSE;
965 3 return TRUE;
966 }
967
968 //=============================================================================
969
970 146 int getRainInterval(int i)
971 //
972 // Input: i = UH group index
973 // Output: returns a time interval (sec)
974 // Purpose: finds rainfall processing time interval for a unit hydrograph group.
975 //
976 {
977 int ri; // rainfal processing time interval for the UH group
978 int tLimb; // duration of a UH's rising & falling limbs
979 int k, m;
980
981 // --- begin with UH group time step equal to wet runoff step
982 146 ri = WetStep;
983
984 // --- examine each UH in the group
985
2/2
✓ Branch 0 taken 1752 times.
✓ Branch 1 taken 146 times.
 1898 for (m=0; m<12; m++)
986 {
987
2/2
✓ Branch 0 taken 5256 times.
✓ Branch 1 taken 1752 times.
 7008 for (k=0; k<3; k++)
988 {
989 // --- make sure the UH exists
990
1/2
✓ Branch 0 taken 5256 times.
✗ Branch 1 not taken.
 5256 if ( UnitHyd[i].tPeak[m][k] > 0 )
991 {
992 // --- reduce time step if rising/falling limb is smaller
993 5256 tLimb = UnitHyd[i].tPeak[m][k];
994 5256 ri = MIN(ri, tLimb);
995 5256 tLimb = UnitHyd[i].tBase[m][k] - tLimb;
996
2/2
✓ Branch 0 taken 5250 times.
✓ Branch 1 taken 6 times.
 5256 if ( tLimb > 0 ) ri = MIN(ri, tLimb);
997 }
998 }
999 }
1000 146 return ri;
1001 }
1002
1003 //=============================================================================
1004
1005 438 int getMaxPeriods(int i, int k)
1006 //
1007 // Input: i = UH group index
1008 // k = UH index
1009 // Output: returns number of past rainfall values
1010 // Purpose: finds number of past rainfall values to save for a UH.
1011 //
1012 {
1013 int m, // month index
1014 n, // number of time periods
1015 nMax, // maximum number of time periods
1016 rainInterval; // rainfall processing interval (sec)
1017
1018 // --- examine each monthly set of UHs
1019 438 rainInterval = UHGroup[i].rainInterval;
1020 438 nMax = 0;
1021
2/2
✓ Branch 0 taken 5256 times.
✓ Branch 1 taken 438 times.
 5694 for (m=0; m<12; m++)
1022 {
1023 // --- compute number of time periods in UH base
1024 5256 n = (UnitHyd[i].tBase[m][k] / rainInterval) + 1;
1025
1026 // --- update number of time periods to be saved
1027 5256 nMax = MAX(n, nMax);
1028 }
1029 438 return nMax;
1030 }
1031
1032 //=============================================================================
1033
1034 3 void initGageData()
1035 //
1036 // Input: none
1037 // Output: none
1038 // Purpose: initializes state of Unit Hydrograph rain gages.
1039 //
1040 {
1041 int i; // unit hyd. index
1042 int g; // rain gage index
1043
1044 // --- first initialize the state of each rain gage
1045
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
 12 for (g=0; g<Nobjects[GAGE]; g++)
1046 {
1047
1/2
✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
 9 if ( Gage[g].tSeries >= 0 )
1048 {
1049 9 table_tseriesInit(&Tseries[Gage[g].tSeries]);
1050 }
1051 9 gage_initState(g);
1052 }
1053
1054 // --- then flag each gage that is used by a Unit Hydrograph set
1055
2/2
✓ Branch 0 taken 146 times.
✓ Branch 1 taken 3 times.
 149 for (i=0; i<Nobjects[UNITHYD]; i++)
1056 {
1057 146 g = UnitHyd[i].rainGage;
1058
1/2
✓ Branch 0 taken 146 times.
✗ Branch 1 not taken.
 146 if ( g >= 0 )
1059 {
1060 // --- if UH's gage uses same time series as a previous gage,
1061 // then assign the latter gage to the UH
1062
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 146 times.
 146 if ( Gage[g].coGage >= 0 )
1063 {
1064 ✗ UnitHyd[i].rainGage = Gage[g].coGage;
1065 }
1066 }
1067 }
1068 3 }
1069
1070 //=============================================================================
1071
1072 3 void initUnitHydData()
1073 //
1074 // Input: none
1075 // Output: none
1076 // Purpose: initializes unit hydrograph processing data.
1077 //
1078 {
1079 int i, // UH group index
1080 j, // node index
1081 k, // UH index
1082 n; // RDII node index
1083 // int g, // rain gage index
1084 int month; // month index
1085
1086 // --- initialize UHGroup entries for each Unit Hydrograph
1087 3 month = datetime_monthOfYear(StartDateTime) - 1;
1088
2/2
✓ Branch 0 taken 146 times.
✓ Branch 1 taken 3 times.
 149 for (i=0; i<Nobjects[UNITHYD]; i++)
1089 {
1090
2/2
✓ Branch 0 taken 438 times.
✓ Branch 1 taken 146 times.
 584 for (k=0; k<3; k++)
1091 {
1092 // --- make the first recorded rainfall begin a new RDII event
1093 // --- (new RDII event occurs when dry period > base of longest UH)
1094 438 UHGroup[i].uh[k].drySeconds =
1095 438 (UHGroup[i].uh[k].maxPeriods * UHGroup[i].rainInterval) + 1;
1096 438 UHGroup[i].uh[k].period = UHGroup[i].uh[k].maxPeriods + 1;
1097 438 UHGroup[i].uh[k].hasPastRain = FALSE;
1098
1099 // --- assign initial abstraction used
1100 438 UHGroup[i].uh[k].iaUsed = UnitHyd[i].iaInit[month][k];
1101 }
1102
1103 // --- initialize gage date to simulation start date
1104 146 UHGroup[i].gageDate = StartDateTime;
1105 146 UHGroup[i].area = 0.0;
1106 146 UHGroup[i].rdii = 0.0;
1107 }
1108
1109 // --- assume each UH group is not used
1110
2/2
✓ Branch 0 taken 146 times.
✓ Branch 1 taken 3 times.
 149 for (i=0; i<Nobjects[UNITHYD]; i++) UHGroup[i].isUsed = FALSE;
1111
1112 // --- look at each node with RDII inflow
1113
2/2
✓ Branch 0 taken 2031 times.
✓ Branch 1 taken 3 times.
 2034 for (n=0; n<NumRdiiNodes; n++)
1114 {
1115 // --- mark as used the UH group associated with the node
1116 2031 j = RdiiNodeIndex[n];
1117 2031 i = Node[j].rdiiInflow->unitHyd;
1118 2031 UHGroup[i].isUsed = TRUE;
1119
1120 // --- add node's sewer area to UH group's area
1121 2031 UHGroup[i].lastDate = StartDateTime;
1122 2031 UHGroup[i].area += Node[j].rdiiInflow->area;
1123 }
1124 3 }
1125
1126 //=============================================================================
1127
1128 3 int openNewRdiiFile()
1129 //
1130 // Input: none
1131 // Output: returns TRUE if successful, FALSE if not
1132 // Purpose: opens a new RDII interface file.
1133 //
1134 {
1135 int j; // node index
1136
1137 // --- create a temporary file name if scratch file being used
1138
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
 3 if ( Frdii.mode == SCRATCH_FILE ) getTempFileName(Frdii.name);
1139
1140 // --- open the RDII file as a formatted text file
1141 3 Frdii.file = fopen(Frdii.name, "w+b");
1142
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if ( Frdii.file == NULL )
1143 {
1144 ✗ return FALSE;
1145 }
1146
1147 // --- write file stamp to RDII file
1148 3 fwrite(FileStamp, sizeof(char), strlen(FileStamp), Frdii.file);
1149
1150 // --- initialize the contents of the file with RDII time step (sec),
1151 // number of RDII nodes, and index of each node
1152 3 fwrite(&RdiiStep, sizeof(INT4), 1, Frdii.file);
1153 3 fwrite(&NumRdiiNodes, sizeof(INT4), 1, Frdii.file);
1154
2/2
✓ Branch 0 taken 2308 times.
✓ Branch 1 taken 3 times.
 2311 for (j=0; j<Nobjects[NODE]; j++)
1155 {
1156
2/2
✓ Branch 0 taken 2031 times.
✓ Branch 1 taken 277 times.
 2308 if ( Node[j].rdiiInflow ) fwrite(&j, sizeof(INT4), 1, Frdii.file);
1157 }
1158 3 return TRUE;
1159 }
1160
1161 //=============================================================================
1162
1163 3039 void getRainfall(DateTime currentDate)
1164 //
1165 // Input: currentDate = current calendar date/time
1166 // Output: none
1167 // Purpose: determines rainfall at current RDII processing date.
1168 //
1169 //
1170 {
1171 int j; // UH group index
1172 int k; // UH index
1173 int g; // rain gage index
1174 int i; // past rainfall index
1175 int month; // month of current date
1176 int rainInterval; // rainfall interval (sec)
1177 double rainDepth; // rainfall depth (inches or mm)
1178 double excessDepth; // excess rainfall depth (inches or mm))
1179 DateTime gageDate; // calendar date for rain gage
1180
1181 // --- examine each UH group
1182 3039 month = datetime_monthOfYear(currentDate) - 1;
1183
2/2
✓ Branch 0 taken 3513 times.
✓ Branch 1 taken 3039 times.
 6552 for (g = 0; g < Nobjects[GAGE]; g++) Gage[g].isCurrent = FALSE;
1184
2/2
✓ Branch 0 taken 17138 times.
✓ Branch 1 taken 3039 times.
 20177 for (j = 0; j < Nobjects[UNITHYD]; j++)
1185 {
1186 // --- repeat until gage's date reaches or exceeds current date
1187 17138 g = UnitHyd[j].rainGage;
1188 17138 rainInterval = UHGroup[j].rainInterval;
1189
2/2
✓ Branch 0 taken 20580 times.
✓ Branch 1 taken 17138 times.
 37718 while ( UHGroup[j].gageDate < currentDate )
1190 {
1191 // --- get rainfall volume over gage's recording interval
1192 // at gage'a current date (in original depth units)
1193 20580 gageDate = UHGroup[j].gageDate;
1194 20580 Adjust.rainFactor = Adjust.rain[datetime_monthOfYear(gageDate)-1];
1195
2/2
✓ Branch 0 taken 3504 times.
✓ Branch 1 taken 17076 times.
 20580 if (!Gage[g].isCurrent)
1196 {
1197 3504 gage_setState(g, gageDate);
1198 3504 Gage[g].isCurrent = TRUE;
1199 }
1200 20580 rainDepth = Gage[g].rainfall * (double)rainInterval / 3600.0;
1201
1202 // --- update amount of total rainfall volume (ft3)
1203 20580 TotalRainVol += rainDepth / UCF(RAINDEPTH) * UHGroup[j].area;
1204
1205 // --- compute rainfall excess for each UH in the group
1206
2/2
✓ Branch 0 taken 61740 times.
✓ Branch 1 taken 20580 times.
 82320 for (k=0; k<3; k++)
1207 {
1208 // --- adjust rainfall volume for any initial abstraction
1209 61740 excessDepth = applyIA(j, k, gageDate, rainInterval, rainDepth);
1210
1211 // --- adjust extent of dry period for the UH
1212 61740 updateDryPeriod(j, k, excessDepth, rainInterval);
1213
1214 // --- add rainfall to list of past values,
1215 // wrapping array index if necessary
1216 61740 i = UHGroup[j].uh[k].period;
1217
2/2
✓ Branch 0 taken 766 times.
✓ Branch 1 taken 60974 times.
 61740 if ( i >= UHGroup[j].uh[k].maxPeriods ) i = 0;
1218 61740 UHGroup[j].uh[k].pastRain[i] = excessDepth;
1219 61740 UHGroup[j].uh[k].pastMonth[i] = (char)month;
1220 61740 UHGroup[j].uh[k].period = i + 1;
1221 }
1222
1223 // --- advance rain date by gage recording interval
1224 20580 UHGroup[j].gageDate = datetime_addSeconds(gageDate, rainInterval);
1225 }
1226 }
1227 3039 }
1228
1229 //=============================================================================
1230
1231 61740 double applyIA(int j, int k, DateTime aDate, double dt, double rainDepth)
1232 //
1233 // Input: j = UH group index
1234 // k = unit hydrograph index
1235 // aDate = current date/time
1236 // dt = time interval (sec)
1237 // rainDepth = unadjusted rain depth (in or mm)
1238 // Output: returns rainfall adjusted for initial abstraction (IA)
1239 // Purpose: adjusts rainfall for any initial abstraction and updates the
1240 // amount of available initial abstraction actually used.
1241 //
1242 {
1243 int m;
1244 double ia, netRainDepth;
1245
1246 // --- determine amount of unused IA
1247 61740 m = datetime_monthOfYear(aDate) - 1;
1248 61740 ia = UnitHyd[j].iaMax[m][k] - UHGroup[j].uh[k].iaUsed;
1249
1/2
✓ Branch 0 taken 61740 times.
✗ Branch 1 not taken.
 61740 ia = MAX(ia, 0.0);
1250
1251 // --- case where there's some rainfall
1252
2/2
✓ Branch 0 taken 8262 times.
✓ Branch 1 taken 53478 times.
 61740 if ( rainDepth > 0.0 )
1253 {
1254 // --- reduce rain depth by unused IA
1255 8262 netRainDepth = rainDepth - ia;
1256
2/2
✓ Branch 0 taken 8244 times.
✓ Branch 1 taken 18 times.
 8262 netRainDepth = MAX(netRainDepth, 0.0);
1257
1258 // --- update amount of IA used up
1259 8262 ia = rainDepth - netRainDepth;
1260 8262 UHGroup[j].uh[k].iaUsed += ia;
1261 }
1262
1263 // --- case where there's no rainfall
1264 else
1265 {
1266 // --- recover a portion of the IA already used
1267 53478 UHGroup[j].uh[k].iaUsed -= dt / 86400. * UnitHyd[j].iaRecov[m][k];
1268
2/2
✓ Branch 0 taken 45828 times.
✓ Branch 1 taken 7650 times.
 53478 UHGroup[j].uh[k].iaUsed = MAX(UHGroup[j].uh[k].iaUsed, 0.0);
1269 53478 netRainDepth = 0.0;
1270 }
1271 61740 return netRainDepth;
1272 }
1273
1274 //=============================================================================
1275
1276 61740 void updateDryPeriod(int j, int k, double rainDepth, int rainInterval)
1277 //
1278 // Input: j = UH group index
1279 // k = unit hydrograph index
1280 // rainDepth = excess rain depth (in or mm)
1281 // rainInterval = rainfall time interval (sec)
1282 // Output: none
1283 // Purpose: adjusts the length of the dry period between rainfall events.
1284 //
1285 {
1286 int i;
1287
1288 // --- if rainfall occurs
1289
2/2
✓ Branch 0 taken 8244 times.
✓ Branch 1 taken 53496 times.
 61740 if ( rainDepth > 0.0 )
1290 {
1291 // --- if previous dry period long enough then begin
1292 // new RDII event with time period index set to 0
1293 8244 if ( UHGroup[j].uh[k].drySeconds >= rainInterval *
1294
2/2
✓ Branch 0 taken 391 times.
✓ Branch 1 taken 7853 times.
 8244 UHGroup[j].uh[k].maxPeriods )
1295 {
1296
2/2
✓ Branch 0 taken 457215 times.
✓ Branch 1 taken 391 times.
 457606 for (i=0; i<UHGroup[j].uh[k].maxPeriods; i++)
1297 {
1298 457215 UHGroup[j].uh[k].pastRain[i] = 0.0;
1299 }
1300 391 UHGroup[j].uh[k].period = 0;
1301 }
1302 8244 UHGroup[j].uh[k].drySeconds = 0;
1303 8244 UHGroup[j].uh[k].hasPastRain = TRUE;
1304 }
1305
1306 // --- if no rainfall, update duration of dry period
1307 else
1308 {
1309 53496 UHGroup[j].uh[k].drySeconds += rainInterval;
1310 53496 if ( UHGroup[j].uh[k].drySeconds >=
1311
2/2
✓ Branch 0 taken 38922 times.
✓ Branch 1 taken 14574 times.
 53496 rainInterval * UHGroup[j].uh[k].maxPeriods )
1312 {
1313 38922 UHGroup[j].uh[k].hasPastRain = FALSE;
1314 }
1315 14574 else UHGroup[j].uh[k].hasPastRain = TRUE;
1316 }
1317 61740 }
1318
1319 //=============================================================================
1320
1321 3039 void getUnitHydRdii(DateTime currentDate)
1322 //
1323 // Input: currentDate = current calendar date/time
1324 // Output: none
1325 // Purpose: computes RDII generated by past rainfall for each UH group.
1326 //
1327 {
1328 int j; // UH group index
1329 int k; // UH index
1330 int rainInterval; // rainfall time interval (sec)
1331
1332 // --- examine each UH group
1333
2/2
✓ Branch 0 taken 17138 times.
✓ Branch 1 taken 3039 times.
 20177 for (j=0; j<Nobjects[UNITHYD]; j++)
1334 {
1335 // --- skip calculation if group not used by any RDII node or if
1336 // current date hasn't reached last date RDII was computed
1337
2/2
✓ Branch 0 taken 3671 times.
✓ Branch 1 taken 13467 times.
 17138 if ( !UHGroup[j].isUsed ) continue;
1338
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13467 times.
 13467 if ( currentDate < UHGroup[j].lastDate ) continue;
1339
1340 // --- update date RDII last computed
1341 13467 UHGroup[j].lastDate = UHGroup[j].gageDate;
1342
1343 // --- perform convolution for each UH in the group
1344 13467 rainInterval = UHGroup[j].rainInterval;
1345 13467 UHGroup[j].rdii = 0.0;
1346
2/2
✓ Branch 0 taken 40401 times.
✓ Branch 1 taken 13467 times.
 53868 for (k=0; k<3; k++)
1347 {
1348
2/2
✓ Branch 0 taken 14785 times.
✓ Branch 1 taken 25616 times.
 40401 if ( UHGroup[j].uh[k].hasPastRain )
1349 {
1350 14785 UHGroup[j].rdii += getUnitHydConvol(j, k, rainInterval);
1351 }
1352 }
1353 }
1354 3039 }
1355
1356 //=============================================================================
1357
1358 14785 double getUnitHydConvol(int j, int k, int rainInterval)
1359 //
1360 // Input: j = UH group index
1361 // k = UH index
1362 // rainInterval = rainfall time interval (sec)
1363 // Output: returns a RDII flow value
1364 // Purpose: computes convolution of Unit Hydrographs with past rainfall.
1365 //
1366 {
1367 int i; // previous rainfall period index
1368 int m; // month of year index
1369 int p; // UH time period index
1370 int pMax; // max. number of periods
1371 double t; // UH time value (sec)
1372 double u; // UH ordinate
1373 double v; // rainfall volume
1374 double rdii; // RDII flow
1375 TUHData* uh; // UH data
1376
1377 // --- initialize RDII, rain period index and UH period index
1378 14785 rdii = 0.0;
1379 14785 uh = &UHGroup[j].uh[k];
1380 14785 i = uh->period - 1;
1381
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14785 times.
 14785 if ( i < 0 ) i = uh->maxPeriods - 1;
1382 14785 pMax = uh->maxPeriods;
1383 14785 p = 1;
1384
1385 // --- evaluate each time period of UH's
1386
2/2
✓ Branch 0 taken 13425260 times.
✓ Branch 1 taken 14785 times.
 13440045 while ( p < pMax )
1387 {
1388 // --- if rain period has rainfall
1389 13425260 v = uh->pastRain[i];
1390 13425260 m = uh->pastMonth[i];
1391
2/2
✓ Branch 0 taken 276422 times.
✓ Branch 1 taken 13148838 times.
 13425260 if ( v > 0.0 )
1392 {
1393 // --- find mid-point time of UH period in seconds
1394 276422 t = ((double)(p) - 0.5) * (double)rainInterval;
1395
1396 // --- convolute rain volume with UH ordinate
1397 276422 u = getUnitHydOrd(j, m, k, t) * UnitHyd[j].r[m][k];
1398 276422 rdii += u * v;
1399 }
1400
1401 // --- move to next UH period & previous rainfall period
1402 13425260 p = p + 1;
1403 13425260 i = i - 1;
1404
2/2
✓ Branch 0 taken 14688 times.
✓ Branch 1 taken 13410572 times.
 13425260 if ( i < 0 ) i = uh->maxPeriods - 1;
1405 }
1406 14785 return rdii;
1407 }
1408
1409 //=============================================================================
1410
1411 276422 double getUnitHydOrd(int h, int m, int k, double t)
1412 //
1413 // Input: h = index of UH group
1414 // m = month index
1415 // k = individual UH index
1416 // t = UH time (sec)
1417 // Output: returns ordinate of a unit hydrograph
1418 // Purpose: gets ordinate of a particular unit hydrograph at specified time.
1419 //
1420 {
1421 double qPeak; // peak flow of unit hydrograph
1422 double f; // fraction of time to/from peak on UH
1423 double t1; // time to peak on UH (sec)
1424 double t2; // time after peak on UH (sec)
1425 double tBase; // base time of UH (sec)
1426
1427 // --- return 0 if past end of UH time base
1428 276422 tBase = UnitHyd[h].tBase[m][k];
1429
2/2
✓ Branch 0 taken 2379 times.
✓ Branch 1 taken 274043 times.
 276422 if ( t >= tBase ) return 0.0;
1430
1431 // --- compute peak value of UH in original rainfall units (in/hr or mm/hr)
1432 274043 qPeak = 2. / tBase * 3600.0;
1433
1434 // --- break UH base into times before & after peak flow
1435 274043 t1 = UnitHyd[h].tPeak[m][k];
1436 274043 t2 = tBase - t1;
1437
1438 // --- find UH flow at time t
1439
2/2
✓ Branch 0 taken 133646 times.
✓ Branch 1 taken 140397 times.
 274043 if ( t <= t1 ) f = t / t1;
1440 140397 else f = 1.0 - (t - t1) / t2;
1441
1/2
✓ Branch 0 taken 274043 times.
✗ Branch 1 not taken.
 274043 return MAX(f, 0.0) * qPeak;
1442 }
1443
1444 //=============================================================================
1445
1446 3039 int getNodeRdii()
1447 //
1448 // Input: none
1449 // Output: returns TRUE if any node has RDII inflow, FALSE if not
1450 // Purpose: computes current RDII inflow at each node.
1451 //
1452 {
1453 3039 int hasRdii = FALSE; // true if any node has some RDII
1454 int i; // UH group index
1455 int j; // node index
1456 int n; // number of nodes w/ RDII
1457 double rdii; // RDII flow (cfs)
1458
1459 // --- examine each node w/ RDII data
1460
2/2
✓ Branch 0 taken 163251 times.
✓ Branch 1 taken 3039 times.
 166290 for (n = 0; n < NumRdiiNodes; n++)
1461 {
1462 // --- identify node's index in project's data base
1463 163251 j = RdiiNodeIndex[n];
1464
1465 // --- apply node's sewer area to UH RDII to get node RDII in CFS
1466 163251 i = Node[j].rdiiInflow->unitHyd;
1467 163251 rdii = UHGroup[i].rdii * Node[j].rdiiInflow->area / UCF(RAINFALL);
1468
2/2
✓ Branch 0 taken 127480 times.
✓ Branch 1 taken 35771 times.
 163251 if ( rdii < ZERO_RDII ) rdii = 0.0;
1469 35771 else hasRdii = TRUE;
1470
1471 // --- update total RDII volume
1472 163251 RdiiNodeFlow[n] = (REAL4)rdii;
1473
2/2
✓ Branch 0 taken 35771 times.
✓ Branch 1 taken 127480 times.
 163251 if ( rdii > 0.0 )
1474 {
1475 35771 TotalRdiiVol += rdii * (double)RdiiStep;
1476 }
1477 }
1478 3039 return hasRdii;
1479 }
1480
1481 //=============================================================================
1482
1483 1130 void saveRdiiFlows(DateTime currentDate)
1484 //
1485 // Input: currentDate = current calendar date/time
1486 // Output: none
1487 // Purpose: saves current set of RDII inflows in current flow units to file.
1488 //
1489 {
1490 1130 fwrite(&currentDate, sizeof(DateTime), 1, Frdii.file);
1491 1130 fwrite(RdiiNodeFlow, sizeof(REAL4), NumRdiiNodes, Frdii.file);
1492 1130 }
1493
1494 //=============================================================================
1495
1496 3 void closeRdiiProcessor()
1497 //
1498 // Input: none
1499 // Output: none
1500 // Purpose: closes RDII processing system.
1501 //
1502 {
1503 // --- write rainfall & RDII totals to report file
1504
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if ( !ErrorCode )
1505 {
1506 3 report_writeRdiiStats(TotalRainVol, TotalRdiiVol);
1507 }
1508
1509 // --- free allocated memory and close RDII file
1510 3 freeRdiiMemory();
1511
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if ( Frdii.file ) fclose(Frdii.file);
1512 3 }
1513
1514 //=============================================================================
1515
1516 3 void freeRdiiMemory()
1517 //
1518 // Input: none
1519 // Output: none
1520 // Purpose: frees memory used for RDII processing.
1521 //
1522 {
1523 int i;
1524 int k;
1525
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if ( UHGroup )
1526 {
1527
2/2
✓ Branch 0 taken 146 times.
✓ Branch 1 taken 3 times.
 149 for (i = 0; i < Nobjects[UNITHYD]; i++)
1528 {
1529
2/2
✓ Branch 0 taken 438 times.
✓ Branch 1 taken 146 times.
 584 for (k=0; k<3; k++)
1530 {
1531
1/2
✓ Branch 0 taken 438 times.
✗ Branch 1 not taken.
 438 FREE(UHGroup[i].uh[k].pastRain);
1532
1/2
✓ Branch 0 taken 438 times.
✗ Branch 1 not taken.
 438 FREE(UHGroup[i].uh[k].pastMonth);
1533 }
1534 }
1535
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 FREE(UHGroup);
1536 }
1537
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 FREE(RdiiNodeIndex);
1538
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 FREE(RdiiNodeFlow);
1539 3 }
1540