GCC Code Coverage Report

Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 84.3% 543 / 0 / 644
Functions: 100.0% 30 / 0 / 30
Branches: 63.2% 309 / 0 / 489
climate.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // climate.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 // Climate related functions.
10 //
11 // Update History
12 // ==============
13 // Build 5.1.007:
14 // - NCDC GHCN climate file format added.
15 // - Monthly adjustments for temperature, evaporation & rainfall added.
16 // Build 5.1.008:
17 // - Monthly adjustments for hyd. conductivity added.
18 // - Time series evaporation rates can now vary within a day.
19 // - Evaporation rates are now properly updated when only flow routing
20 // is being simulated.
21 // Build 5.1.010:
22 // - Hargreaves evaporation now computed using 7-day average temperatures.
23 // Build 5.1.011:
24 // - Monthly adjustment for hyd. conductivity <= 0 is ignored.
25 // Build 5.1.013:
26 // - Reads names of monthly adjustment patterns for various parameters
27 // of a subcatchment from the [ADJUSTMENTS] section of input file.
28 // Build 5.2.0:
29 // - Reads temperature units for use with GHCND climate files.
30 // - Support added for relative file names.
31 ///-----------------------------------------------------------------------------
32 #define _CRT_SECURE_NO_DEPRECATE
33
34 #include <stdlib.h>
35 #include <stdio.h>
36 #include <string.h>
37 #include <math.h>
38 #include "headers.h"
39
40 //-----------------------------------------------------------------------------
41 // Constants
42 //-----------------------------------------------------------------------------
43 enum ClimateFileFormats {UNKNOWN_FORMAT,
44 USER_PREPARED, // SWMM 5's own user format
45 GHCND, // NCDC GHCN Daily format
46 TD3200, // NCDC TD3200 format
47 DLY0204}; // Canadian DLY02 or DLY04 format
48 static const int MAXCLIMATEVARS = 4;
49 static const int MAXDAYSPERMONTH = 32;
50
51 // These variables are used when processing climate files.
52 enum ClimateVarType {TMIN, TMAX, EVAP, WIND};
53 enum WindSpeedType {WDMV, AWND};
54 enum TempUnitsType {DEG_C10, DEG_C, DEG_F};
55 static char* ClimateVarWords[] = {"TMIN", "TMAX", "EVAP", "WDMV", "AWND",
56 NULL};
57 static char* TempUnitsWords[] = {"C10", "C", "F", NULL};
58
59 //-----------------------------------------------------------------------------
60 // Data Structures
61 //-----------------------------------------------------------------------------
62 typedef struct
63 {
64 double tAve; // moving avg. for daily temperature (deg F)
65 double tRng; // moving avg. for daily temp. range (deg F)
66 double ta[7]; // data window for tAve
67 double tr[7]; // data window for tRng
68 int count; // length of moving average window
69 int maxCount; // maximum length of moving average window
70 int front; // index of front of moving average window
71 } TMovAve;
72
73
74 //-----------------------------------------------------------------------------
75 // Shared variables
76 //-----------------------------------------------------------------------------
77 // Temperature variables
78 static double Tmin; // min. daily temperature (deg F)
79 static double Tmax; // max. daily temperature (deg F)
80 static double Trng; // 1/2 range of daily temperatures
81 static double Trng1; // prev. max - current min. temp.
82 static double Tave; // average daily temperature (deg F)
83 static double Hrsr; // time of min. temp. (hrs)
84 static double Hrss; // time of max. temp (hrs)
85 static double Hrday; // avg. of min/max temp times
86 static double Dhrdy; // hrs. between min. & max. temp. times
87 static double Dydif; // hrs. between max. & min. temp. times
88 static DateTime LastDay; // date of last day with temp. data
89 static TMovAve Tma; // moving average of daily temperatures
90
91 // Evaporation variables
92 static DateTime NextEvapDate; // next date when evap. rate changes
93 static double NextEvapRate; // next evaporation rate (user units)
94
95 // Climate file variables
96 static int FileFormat; // file format (see ClimateFileFormats)
97 static int FileYear; // current year of file data
98 static int FileMonth; // current month of year of file data
99 static int FileDay; // current day of month of file data
100 static int FileLastDay; // last day of current month of file data
101 static int FileElapsedDays; // number of days read from file
102 static double FileValue[4]; // current day's values of climate data
103 static double FileData[4][32]; // month's worth of daily climate data
104 static char FileLine[MAXLINE+1]; // line from climate data file
105
106 static int FileFieldPos[4]; // start of data fields for file record
107 static int FileDateFieldPos; // start of date field for file record
108 static int FileWindType; // wind speed type
109 static int FileTempUnits; // GHCND file temperature units (C10, C or F)
110
111 //-----------------------------------------------------------------------------
112 // External functions (defined in funcs.h)
113 //-----------------------------------------------------------------------------
114 // climate_readParams // called by input_parseLine
115 // climate_readEvapParams // called by input_parseLine
116 // climate_validate // called by project_validate
117 // climate_openFile // called by runoff_open
118 // climate_initState // called by project_init
119 // climate_setState // called by runoff_execute
120 // climate_getNextEvapDate // called by runoff_getTimeStep
121
122 //-----------------------------------------------------------------------------
123 // Local functions
124 //-----------------------------------------------------------------------------
125 static int getFileFormat(void);
126 static void readFileLine(int *year, int *month);
127 static void readUserFileLine(int *year, int *month);
128 static void readTD3200FileLine(int *year, int *month);
129 static void readDLY0204FileLine(int *year, int *month);
130 static void readFileValues(void);
131
132 static void setNextEvapDate(DateTime thedate);
133 static void setEvap(DateTime theDate);
134 static void setTemp(DateTime theDate);
135 static void setWind(DateTime theDate);
136 static void updateTempTimes(int day);
137 static void updateTempMoveAve(double tmin, double tmax);
138 static double getTempEvap(int day, double ta, double tr);
139
140 static void updateFileValues(DateTime theDate);
141 static void parseUserFileLine(void);
142 static void parseTD3200FileLine(void);
143 static void parseDLY0204FileLine(void);
144 static void setTD3200FileValues(int param);
145
146 static int isGhcndFormat(char* line);
147 static void readGhcndFileLine(int *year, int *month);
148 static void parseGhcndFileLine(void);
149 static double convertGhcndValue(int var, double v);
150
151 //=============================================================================
152
153 48 int climate_readParams(char* tok[], int ntoks)
154 //
155 // Input: tok[] = array of string tokens
156 // ntoks = number of tokens
157 // Output: returns error code
158 // Purpose: reads climate/temperature parameters from input line of data
159 //
160 // Format of data can be
161 // TIMESERIES name
162 // FILE name (start) (units)
163 // WINDSPEED MONTHLY v1 v2 ... v12
164 // WINDSPEED FILE
165 // SNOWMELT v1 v2 ... v6
166 // ADC IMPERV/PERV v1 v2 ... v10
167 //
168 {
169 int i, j, k;
170 double x[6], y;
171 DateTime aDate;
172 char fname[MAXFNAME + 1];
173
174 // --- identify keyword
175 48 k = findmatch(tok[0], TempKeyWords);
176
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
 48 if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[0]);
177
5/6
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 48 switch (k)
178 {
179 1 case 0: // Time series name
180 // --- check that time series name exists
181
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( ntoks < 2 ) return error_setInpError(ERR_ITEMS, "");
182 1 i = project_findObject(TSERIES, tok[1]);
183
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( i < 0 ) return error_setInpError(ERR_NAME, tok[1]);
184
185 // --- record the time series as being the data source for temperature
186 1 Temp.dataSource = TSERIES_TEMP;
187 1 Temp.tSeries = i;
188 1 Tseries[i].refersTo = TSERIES_TEMP;
189 1 break;
190
191 11 case 1: // Climate file
192 // --- record file as being source of temperature data
193
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if ( ntoks < 2 ) return error_setInpError(ERR_ITEMS, "");
194 11 Temp.dataSource = FILE_TEMP;
195
196 // --- save name and usage mode of external climate file
197 11 Fclimate.mode = USE_FILE;
198 11 sstrncpy(fname, tok[1], MAXFNAME);
199 11 sstrncpy(Fclimate.name, addAbsolutePath(fname), MAXFNAME);
200
201 // --- save starting date to read from file if one is provided
202 11 Temp.fileStartDate = NO_DATE;
203
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 7 times.
 11 if ( ntoks > 2 )
204 {
205
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if ( *tok[2] != '*')
206 {
207
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( !datetime_strToDate(tok[2], &aDate) )
208 ✗ return error_setInpError(ERR_DATETIME, tok[2]);
209 4 Temp.fileStartDate = aDate;
210 }
211 }
212
213 // --- file temperature units
214 11 FileTempUnits = DEG_F;
215
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if (UnitSystem == SI)
216 ✗ FileTempUnits = DEG_C;
217
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 7 times.
 11 if (ntoks > 3)
218 {
219 4 i = findmatch(tok[3], TempUnitsWords);
220
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (i < 0)
221 ✗ return error_setInpError(ERR_KEYWORD, tok[3]);
222 4 FileTempUnits = i;
223 }
224 11 break;
225
226 11 case 2: // Wind speeds
227 // --- check if wind speeds will be supplied from climate file
228
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 10 times.
 11 if ( strcomp(tok[1], w_FILE) )
229 {
230 1 Wind.type = FILE_WIND;
231 }
232
233 // --- otherwise read 12 monthly avg. wind speed values
234 else
235 {
236
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( ntoks < 14 ) return error_setInpError(ERR_ITEMS, "");
237 10 Wind.type = MONTHLY_WIND;
238
2/2
✓ Branch 0 taken 120 times.
✓ Branch 1 taken 10 times.
 130 for (i=0; i<12; i++)
239 {
240
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 120 times.
 120 if ( !getDouble(tok[i+2], &y) )
241 ✗ return error_setInpError(ERR_NUMBER, tok[i+2]);
242 120 Wind.aws[i] = y;
243 }
244 }
245 11 break;
246
247 9 case 3: // Snowmelt params
248
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 if ( ntoks < 7 ) return error_setInpError(ERR_ITEMS, "");
249
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 9 times.
 63 for (i=1; i<7; i++)
250 {
251
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 54 times.
 54 if ( !getDouble(tok[i], &x[i-1]) )
252 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
253 }
254 // --- convert deg. C to deg. F for snowfall temperature
255
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 if ( UnitSystem == SI ) x[0] = 9./5.*x[0] + 32.0;
256 9 Snow.snotmp = x[0];
257 9 Snow.tipm = x[1];
258 9 Snow.rnm = x[2];
259 9 Temp.elev = x[3] / UCF(LENGTH);
260 9 Temp.anglat = x[4];
261 9 Temp.dtlong = x[5] / 60.0;
262 9 break;
263
264 16 case 4: // Areal Depletion Curve data
265 // --- check if data is for impervious or pervious areas
266
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
 16 if ( ntoks < 12 ) return error_setInpError(ERR_ITEMS, "");
267
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 8 times.
 16 if ( match(tok[1], w_IMPERV) ) i = 0;
268
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 else if ( match(tok[1], w_PERV) ) i = 1;
269 ✗ else return error_setInpError(ERR_KEYWORD, tok[1]);
270
271 // --- read 10 fractional values
272
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 16 times.
 176 for (j=0; j<10; j++)
273 {
274
3/6
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 160 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 160 times.
 160 if ( !getDouble(tok[j+2], &y) || y < 0.0 || y > 1.0 )
275 ✗ return error_setInpError(ERR_NUMBER, tok[j+2]);
276 160 Snow.adc[i][j] = y;
277 }
278 16 break;
279 }
280 48 return 0;
281 }
282
283 //=============================================================================
284
285 92 int climate_readEvapParams(char* tok[], int ntoks)
286 //
287 // Input: tok[] = array of string tokens
288 // ntoks = number of tokens
289 // Output: returns error code
290 // Purpose: reads evaporation parameters from input line of data.
291 //
292 // Data formats are:
293 // CONSTANT value
294 // MONTHLY v1 ... v12
295 // TIMESERIES name
296 // TEMPERATURE
297 // FILE (v1 ... v12)
298 // RECOVERY name
299 // DRY_ONLY YES/NO
300 //
301 {
302 int i, k;
303 double x;
304
305 // --- find keyword indicating what form the evaporation data is in
306 92 k = findmatch(tok[0], EvapTypeWords);
307
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
 92 if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[0]);
308
309 // --- check for RECOVERY pattern data
310
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 90 times.
 92 if ( k == RECOVERY )
311 {
312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 2 ) return error_setInpError(ERR_ITEMS, "");
313 2 i = project_findObject(TIMEPATTERN, tok[1]);
314
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( i < 0 ) return error_setInpError(ERR_NAME, tok[1]);
315 2 Evap.recoveryPattern = i;
316 2 return 0;
317 }
318
319 // --- check for no evaporation in wet periods
320
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 46 times.
 90 if ( k == DRYONLY )
321 {
322
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
 44 if ( ntoks < 2 ) return error_setInpError(ERR_ITEMS, "");
323
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 if ( strcomp(tok[1], w_NO ) ) Evap.dryOnly = FALSE;
324 ✗ else if ( strcomp(tok[1], w_YES ) ) Evap.dryOnly = TRUE;
325 ✗ else return error_setInpError(ERR_KEYWORD, tok[1]);
326 44 return 0;
327 }
328
329 // --- process data depending on its form
330 46 Evap.type = k;
331
3/4
✓ Branch 0 taken 37 times.
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 37 times.
 46 if ( k != TEMPERATURE_EVAP && ntoks < 2 )
332 ✗ return error_setInpError(ERR_ITEMS, "");
333
5/5
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 9 times.
 46 switch ( k )
334 {
335 32 case CONSTANT_EVAP:
336 // --- for constant evap., fill monthly avg. values with same number
337
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 32 times.
 32 if ( !getDouble(tok[1], &x) )
338 ✗ return error_setInpError(ERR_NUMBER, tok[1]);
339
2/2
✓ Branch 0 taken 384 times.
✓ Branch 1 taken 32 times.
 416 for (i=0; i<12; i++) Evap.monthlyEvap[i] = x;
340 32 break;
341
342 2 case MONTHLY_EVAP:
343 // --- for monthly evap., read a value for each month of year
344
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 13 ) return error_setInpError(ERR_ITEMS, "");
345
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 2 times.
 26 for ( i=0; i<12; i++)
346
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if ( !getDouble(tok[i+1], &Evap.monthlyEvap[i]) )
347 ✗ return error_setInpError(ERR_NUMBER, tok[i+1]);
348 2 break;
349
350 1 case TIMESERIES_EVAP:
351 // --- for time series evap., read name of time series
352 1 i = project_findObject(TSERIES, tok[1]);
353
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( i < 0 ) return error_setInpError(ERR_NAME, tok[1]);
354 1 Evap.tSeries = i;
355 1 Tseries[i].refersTo = TIMESERIES_EVAP;
356 1 break;
357
358 2 case FILE_EVAP:
359 // --- for evap. from climate file, read monthly pan coeffs.
360 // if they are provided (default values are 1.0)
361
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( ntoks > 1 )
362 {
363
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 13 ) return error_setInpError(ERR_ITEMS, "");
364
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 2 times.
 26 for (i=0; i<12; i++)
365 {
366
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if ( !getDouble(tok[i+1], &Evap.panCoeff[i]) )
367 ✗ return error_setInpError(ERR_NUMBER, tok[i+1]);
368 }
369 }
370 2 break;
371 }
372 46 return 0;
373 }
374
375 //=============================================================================
376
377 8 int climate_readAdjustments(char* tok[], int ntoks)
378 //
379 // Input: tok[] = array of string tokens
380 // ntoks = number of tokens
381 // Output: returns error code
382 // Purpose: reads adjustments to monthly evaporation or rainfall
383 // from input line of data.
384 //
385 // Data formats are:
386 // TEMPERATURE v1 ... v12
387 // EVAPORATION v1 ... v12
388 // RAINFALL v1 ... v12
389 // CONDUCTIVITY v1 ... v12
390 // N-PERV subcatchID patternID
391 // DSTORE subcatchID patternID
392 // INFIL subcatchID patternID
393 {
394 int i, j;
395
396
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if (ntoks == 1) return 0;
397
398
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 6 times.
 8 if ( match(tok[0], "TEMP") )
399 {
400
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 13 ) return error_setInpError(ERR_ITEMS, "");
401
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 2 times.
 26 for (i = 1; i < 13; i++)
402 {
403
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if ( !getDouble(tok[i], &Adjust.temp[i-1]) )
404 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
405 }
406 2 return 0;
407 }
408
409
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
 6 if ( match(tok[0], "EVAP") )
410 {
411
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 13 ) return error_setInpError(ERR_ITEMS, "");
412
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 2 times.
 26 for (i = 1; i < 13; i++)
413 {
414
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if ( !getDouble(tok[i], &Adjust.evap[i-1]) )
415 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
416 }
417 2 return 0;
418 }
419
420
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 if ( match(tok[0], "RAIN") )
421 {
422
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 13 ) return error_setInpError(ERR_ITEMS, "");
423
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 2 times.
 26 for (i = 1; i < 13; i++)
424 {
425
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if ( !getDouble(tok[i], &Adjust.rain[i-1]) )
426 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
427 }
428 2 return 0;
429 }
430
431
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 if ( match(tok[0], "CONDUCT") )
432 {
433
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 13 ) return error_setInpError(ERR_ITEMS, "");
434
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 2 times.
 26 for (i = 1; i < 13; i++)
435 {
436
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if ( !getDouble(tok[i], &Adjust.hydcon[i-1]) )
437 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
438
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
 24 if ( Adjust.hydcon[i-1] <= 0.0 ) Adjust.hydcon[i-1] = 1.0;
439 }
440 2 return 0;
441 }
442
443 ✗ if ( match(tok[0], "N-PERV") )
444 {
445 ✗ if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
446 ✗ i = project_findObject(SUBCATCH, tok[1]);
447 ✗ if (i < 0) return error_setInpError(ERR_NAME, tok[1]);
448 ✗ j = project_findObject(TIMEPATTERN, tok[2]);
449 ✗ if (j < 0) return error_setInpError(ERR_NAME, tok[2]);
450 ✗ Subcatch[i].nPervPattern = j;
451 ✗ return 0;
452 }
453
454 ✗ if ( match(tok[0], "DSTORE") )
455 {
456 ✗ if (ntoks < 3) return error_setInpError(ERR_ITEMS, "");
457 ✗ i = project_findObject(SUBCATCH, tok[1]);
458 ✗ if (i < 0) return error_setInpError(ERR_NAME, tok[1]);
459 ✗ j = project_findObject(TIMEPATTERN, tok[2]);
460 ✗ if (j < 0) return error_setInpError(ERR_NAME, tok[2]);
461 ✗ Subcatch[i].dStorePattern = j;
462 ✗ return 0;
463 }
464
465 ✗ if (match(tok[0], "INFIL"))
466 {
467 ✗ if (ntoks < 3) return error_setInpError(ERR_ITEMS, "");
468 ✗ i = project_findObject(SUBCATCH, tok[1]);
469 ✗ if (i < 0) return error_setInpError(ERR_NAME, tok[1]);
470 ✗ j = project_findObject(TIMEPATTERN, tok[2]);
471 ✗ if (j < 0) return error_setInpError(ERR_NAME, tok[2]);
472 ✗ Subcatch[i].infilPattern = j;
473 ✗ return 0;
474 }
475 ✗ return error_setInpError(ERR_KEYWORD, tok[0]);
476 }
477
478 //=============================================================================
479
480 58 void climate_validate()
481 //
482 // Input: none
483 // Output: none
484 // Purpose: validates climatological variables
485 //
486 {
487 int i;
488 double a, z, pa;
489
490 // --- check if climate data comes from external data file
491
4/4
✓ Branch 0 taken 57 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 56 times.
✓ Branch 3 taken 1 time.
 58 if ( Wind.type == FILE_WIND || Evap.type == FILE_EVAP ||
492
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 47 times.
 56 Evap.type == TEMPERATURE_EVAP )
493 {
494
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if ( Fclimate.mode == NO_FILE )
495 {
496 ✗ report_writeErrorMsg(ERR_NO_CLIMATE_FILE, "");
497 }
498 }
499
500 // --- open the climate data file
501
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 47 times.
 58 if ( Fclimate.mode == USE_FILE ) climate_openFile();
502
503 // --- snow melt parameters tipm & rnm must be fractions
504
1/2
✓ Branch 0 taken 58 times.
✗ Branch 1 not taken.
 58 if ( Snow.tipm < 0.0 ||
505
1/2
✓ Branch 0 taken 58 times.
✗ Branch 1 not taken.
 58 Snow.tipm > 1.0 ||
506
1/2
✓ Branch 0 taken 58 times.
✗ Branch 1 not taken.
 58 Snow.rnm < 0.0 ||
507
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 58 times.
 58 Snow.rnm > 1.0 ) report_writeErrorMsg(ERR_SNOWMELT_PARAMS, "");
508
509 // --- latitude should be between -90 & 90 degrees
510 58 a = Temp.anglat;
511
2/4
✓ Branch 0 taken 58 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 58 times.
 58 if ( a <= -89.99 ||
512 ✗ a >= 89.99 ) report_writeErrorMsg(ERR_SNOWMELT_PARAMS, "");
513 58 else Temp.tanAnglat = tan(a * PI / 180.0);
514
515 // --- compute psychrometric constant
516 58 z = Temp.elev / 1000.0;
517
2/2
✓ Branch 0 taken 49 times.
✓ Branch 1 taken 9 times.
 58 if ( z <= 0.0 ) pa = 29.9;
518 9 else pa = 29.9 - 1.02*z + 0.0032*pow(z, 2.4); // atmos. pressure
519 58 Temp.gamma = 0.000359 * pa;
520
521 // --- convert units of monthly temperature & evap adjustments
522
2/2
✓ Branch 0 taken 696 times.
✓ Branch 1 taken 58 times.
 754 for (i = 0; i < 12; i++)
523 {
524
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 660 times.
 696 if (UnitSystem == SI) Adjust.temp[i] *= 9.0/5.0;
525 696 Adjust.evap[i] /= UCF(EVAPRATE);
526 }
527 58 }
528
529 //=============================================================================
530
531 11 void climate_openFile()
532 //
533 // Input: none
534 // Output: none
535 // Purpose: opens a climate file and reads in first set of values.
536 //
537 {
538 int i, m, y;
539
540 // --- open the file
541
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 11 times.
 11 if ( (Fclimate.file = fopen(Fclimate.name, "rt")) == NULL )
542 {
543 ✗ report_writeErrorMsg(ERR_CLIMATE_FILE_OPEN, Fclimate.name);
544 ✗ return;
545 }
546
547 // --- initialize values of file's climate variables
548 // (Temp.ta was previously initialized in project.c)
549 11 FileValue[TMIN] = Temp.ta;
550 11 FileValue[TMAX] = Temp.ta;
551 11 FileValue[EVAP] = 0.0;
552 11 FileValue[WIND] = 0.0;
553
554 // --- find climate file's format
555 11 FileFormat = getFileFormat();
556
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if ( FileFormat == UNKNOWN_FORMAT )
557 {
558 ✗ report_writeErrorMsg(ERR_CLIMATE_FILE_READ, Fclimate.name);
559 ✗ return;
560 }
561
562 // --- position file to begin reading climate file at either user-specified
563 // month/year or at start of simulation period.
564 11 rewind(Fclimate.file);
565 11 sstrncpy(FileLine, "", 0);
566
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 4 times.
 11 if ( Temp.fileStartDate == NO_DATE )
567 7 datetime_decodeDate(StartDate, &FileYear, &FileMonth, &FileDay);
568 else
569 4 datetime_decodeDate(Temp.fileStartDate, &FileYear, &FileMonth, &FileDay);
570
1/2
✓ Branch 1 taken 29367 times.
✗ Branch 2 not taken.
 29367 while ( !feof(Fclimate.file) )
571 {
572 29367 sstrncpy(FileLine, "", 0);
573 29367 readFileLine(&y, &m);
574
4/4
✓ Branch 0 taken 346 times.
✓ Branch 1 taken 29021 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 335 times.
 29367 if ( y == FileYear && m == FileMonth ) break;
575 }
576
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 11 times.
 11 if ( feof(Fclimate.file) )
577 {
578 ✗ report_writeErrorMsg(ERR_CLIMATE_END_OF_FILE, Fclimate.name);
579 ✗ return;
580 }
581
582 // --- initialize file dates and current climate variable values
583
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if ( !ErrorCode )
584 {
585 11 FileElapsedDays = 0;
586 11 FileLastDay = datetime_daysPerMonth(FileYear, FileMonth);
587 11 readFileValues();
588
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 11 times.
 55 for (i=TMIN; i<=WIND; i++)
589 {
590
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 24 times.
 44 if ( FileData[i][FileDay] == MISSING ) continue;
591 24 FileValue[i] = FileData[i][FileDay];
592 }
593 }
594 }
595
596 //=============================================================================
597
598 58 void climate_initState()
599 //
600 // Input: none
601 // Output: none
602 // Purpose: initializes climate state variables.
603 //
604 {
605 58 LastDay = NO_DATE;
606 58 Temp.tmax = MISSING;
607 58 Snow.removed = 0.0;
608 58 NextEvapDate = StartDate;
609 58 NextEvapRate = 0.0;
610
611 // --- initialize variables for time series evaporation
612
3/4
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 57 times.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 58 if ( Evap.type == TIMESERIES_EVAP && Evap.tSeries >= 0 )
613 {
614 // --- initialize NextEvapDate & NextEvapRate to first entry of
615 // time series whose date <= the simulation start date
616 1 table_getFirstEntry(&Tseries[Evap.tSeries],
617 &NextEvapDate, &NextEvapRate);
618
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( NextEvapDate < StartDate )
619 {
620 ✗ setNextEvapDate(StartDate);
621 }
622 1 Evap.rate = NextEvapRate / UCF(EVAPRATE);
623
624 // --- find the next time evaporation rates change after this
625 1 setNextEvapDate(NextEvapDate);
626 }
627
628 // --- initialize variables for temperature evaporation
629
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 49 times.
 58 if ( Evap.type == TEMPERATURE_EVAP )
630 {
631 9 Tma.maxCount = sizeof(Tma.ta) / sizeof(double);
632 9 Tma.count = 0;
633 9 Tma.front = 0;
634 9 Tma.tAve = 0.0;
635 9 Tma.tRng = 0.0;
636 }
637 58 }
638
639 //=============================================================================
640
641 686476 void climate_setState(DateTime theDate)
642 //
643 // Input: theDate = simulation date
644 // Output: none
645 // Purpose: sets climate variables for current date.
646 //
647 {
648
2/2
✓ Branch 0 taken 11592 times.
✓ Branch 1 taken 674884 times.
 686476 if ( Fclimate.mode == USE_FILE ) updateFileValues(theDate);
649
2/2
✓ Branch 0 taken 12012 times.
✓ Branch 1 taken 674464 times.
 686476 if ( Temp.dataSource != NO_TEMP ) setTemp(theDate);
650 686476 setEvap(theDate);
651 686476 setWind(theDate);
652 686476 Adjust.rainFactor = Adjust.rain[datetime_monthOfYear(theDate)-1];
653 686476 Adjust.hydconFactor = Adjust.hydcon[datetime_monthOfYear(theDate)-1];
654 686476 setNextEvapDate(theDate);
655 686476 }
656
657 //=============================================================================
658
659 33090 DateTime climate_getNextEvapDate()
660 //
661 // Input: none
662 // Output: returns the current value of NextEvapDate
663 // Purpose: gets the next date when evaporation rate changes.
664 //
665 {
666 33090 return NextEvapDate;
667 }
668
669 //=============================================================================
670
671 686477 void setNextEvapDate(DateTime theDate)
672 //
673 // Input: theDate = current simulation date
674 // Output: sets a new value for NextEvapDate
675 // Purpose: finds date for next change in evaporation after the current date.
676 //
677 {
678 int yr, mon, day, k;
679 double d, e;
680
681 // --- do nothing if current date hasn't reached the current next date
682
2/2
✓ Branch 0 taken 686411 times.
✓ Branch 1 taken 66 times.
 686477 if ( NextEvapDate > theDate ) return;
683
684
5/5
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 9 times.
 66 switch ( Evap.type )
685 {
686 // --- for constant evaporation, use a next date far in the future
687 44 case CONSTANT_EVAP:
688 44 NextEvapDate = theDate + 365.;
689 44 break;
690
691 // --- for monthly evaporation, use the start of the next month
692 2 case MONTHLY_EVAP:
693 2 datetime_decodeDate(theDate, &yr, &mon, &day);
694
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( mon == 12 )
695 {
696 ✗ mon = 1;
697 ✗ yr++;
698 }
699 2 else mon++;
700 2 NextEvapDate = datetime_encodeDate(yr, mon, 1);
701 2 break;
702
703 // --- for time series evaporation, find the next entry in the
704 // series on or after the current date
705 3 case TIMESERIES_EVAP:
706 3 k = Evap.tSeries;
707
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if ( k >= 0 )
708 {
709 3 NextEvapDate = theDate + 365.;
710
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 time.
 3 while ( table_getNextEntry(&Tseries[k], &d, &e) &&
711
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 d <= EndDateTime )
712 {
713
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( d >= theDate )
714 {
715 2 NextEvapDate = d;
716 2 NextEvapRate = e;
717 2 break;
718 }
719 }
720 }
721 3 break;
722
723 // --- for climate file daily evaporation, use the next day
724 8 case FILE_EVAP:
725 8 NextEvapDate = floor(theDate) + 1.0;
726 8 break;
727
728 9 default: NextEvapDate = theDate + 365.;
729 }
730 }
731
732 //=============================================================================
733
734 11592 void updateFileValues(DateTime theDate)
735 //
736 // Input: theDate = current simulation date
737 // Output: none
738 // Purpose: updates daily climate variables for new day or reads in
739 // another month worth of values if a new month begins.
740 //
741 // NOTE: counters FileElapsedDays, FileDay, FileMonth, FileYear and
742 // FileLastDay were initialized in climate_openFile().
743 //
744 {
745 int i;
746 int deltaDays;
747
748 // --- see if a new day has begun
749 11592 deltaDays = (int)(floor(theDate) - floor(StartDateTime));
750
2/2
✓ Branch 0 taken 41 times.
✓ Branch 1 taken 11551 times.
 11592 if ( deltaDays > FileElapsedDays )
751 {
752 // --- advance day counters
753 41 FileElapsedDays++;
754 41 FileDay++;
755
756 // --- see if new month of data needs to be read from file
757
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 40 times.
 41 if ( FileDay > FileLastDay )
758 {
759 1 FileMonth++;
760
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if ( FileMonth > 12 )
761 {
762 1 FileMonth = 1;
763 1 FileYear++;
764 }
765 1 readFileValues();
766 1 FileDay = 1;
767 1 FileLastDay = datetime_daysPerMonth(FileYear, FileMonth);
768 }
769
770 // --- set climate variables for new day
771
2/2
✓ Branch 0 taken 164 times.
✓ Branch 1 taken 41 times.
 205 for (i=TMIN; i<=WIND; i++)
772 {
773 // --- no change in current value if its missing
774
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 88 times.
 164 if ( FileData[i][FileDay] == MISSING ) continue;
775 88 FileValue[i] = FileData[i][FileDay];
776 }
777 }
778 11592 }
779
780 //=============================================================================
781
782 12012 void setTemp(DateTime theDate)
783 //
784 // Input: theDate = simulation date
785 // Output: none
786 // Purpose: updates temperatures for new simulation date.
787 //
788 {
789 int j; // snow data object index
790 int k; // time series index
791 int mon; // month of year
792 int day; // day of year
793 DateTime theDay; // calendar day
794 double hour; // hour of day
795 double tmp; // temporary temperature
796
797 // --- see if a new day has started
798 12012 mon = datetime_monthOfYear(theDate);
799 12012 theDay = floor(theDate);
800
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 11958 times.
 12012 if ( theDay > LastDay )
801 {
802 // --- update min. & max. temps & their time of day
803 54 day = datetime_dayOfYear(theDate);
804
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 2 times.
 54 if ( Temp.dataSource == FILE_TEMP )
805 {
806 52 Tmin = FileValue[TMIN] + Adjust.temp[mon-1];
807 52 Tmax = FileValue[TMAX] + Adjust.temp[mon-1];
808
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if ( Tmin > Tmax )
809 {
810 ✗ tmp = Tmin;
811 ✗ Tmin = Tmax;
812 ✗ Tmax = tmp;
813 }
814 52 updateTempTimes(day);
815
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 8 times.
 52 if ( Evap.type == TEMPERATURE_EVAP )
816 {
817 44 updateTempMoveAve(Tmin, Tmax);
818 44 FileValue[EVAP] = getTempEvap(day, Tma.tAve, Tma.tRng);
819 }
820 }
821
822 // --- compute snow melt coefficients based on day of year
823 54 Snow.season = sin(0.0172615*(day-81.0));
824
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 54 times.
 94 for (j=0; j<Nobjects[SNOWMELT]; j++)
825 {
826 40 snow_setMeltCoeffs(j, Snow.season);
827 }
828
829 // --- update date of last day analyzed
830 54 LastDay = theDate;
831 }
832
833 // --- for min/max daily temps. from climate file,
834 // compute hourly temp. by sinusoidal interp.
835
2/2
✓ Branch 0 taken 11592 times.
✓ Branch 1 taken 420 times.
 12012 if ( Temp.dataSource == FILE_TEMP )
836 {
837 11592 hour = (theDate - theDay) * 24.0;
838
2/2
✓ Branch 0 taken 3183 times.
✓ Branch 1 taken 8409 times.
 11592 if ( hour < Hrsr )
839 3183 Temp.ta = Tmin + Trng1/2.0 * sin(PI/Dydif * (Hrsr - hour));
840
3/4
✓ Branch 0 taken 8409 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3429 times.
✓ Branch 3 taken 4980 times.
 8409 else if ( hour >= Hrsr && hour <= Hrss )
841 3429 Temp.ta = Tave + Trng * sin(PI/Dhrdy * (Hrday - hour));
842 else
843 4980 Temp.ta = Tmax - Trng * sin(PI/Dydif * (hour - Hrss));
844 }
845
846 // --- for user-supplied temperature time series,
847 // get temperature value from time series
848
2/2
✓ Branch 0 taken 420 times.
✓ Branch 1 taken 11592 times.
 12012 if ( Temp.dataSource == TSERIES_TEMP )
849 {
850 420 k = Temp.tSeries;
851
1/2
✓ Branch 0 taken 420 times.
✗ Branch 1 not taken.
 420 if ( k >= 0)
852 {
853 420 Temp.ta = table_tseriesLookup(&Tseries[k], theDate, TRUE);
854
855 // --- convert from deg. C to deg. F if need be
856
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 420 times.
 420 if ( UnitSystem == SI )
857 {
858 ✗ Temp.ta = (9./5.) * Temp.ta + 32.0;
859 }
860
861 // --- apply climate change adjustment factor
862 420 Temp.ta += Adjust.temp[mon-1];
863 }
864 }
865
866 // --- compute saturation vapor pressure
867 12012 Temp.ea = 8.1175e6 * exp(-7701.544 / (Temp.ta + 405.0265) );
868 12012 }
869
870 //=============================================================================
871
872 686476 void setEvap(DateTime theDate)
873 //
874 // Input: theDate = simulation date
875 // Output: none
876 // Purpose: sets evaporation rate (ft/sec) for a specified date.
877 //
878 {
879 int k;
880 686476 int mon = datetime_monthOfYear(theDate);
881
882
5/6
✓ Branch 0 taken 674218 times.
✓ Branch 1 taken 90 times.
✓ Branch 2 taken 576 times.
✓ Branch 3 taken 833 times.
✓ Branch 4 taken 10759 times.
✗ Branch 5 not taken.
 686476 switch ( Evap.type )
883 {
884 674218 case CONSTANT_EVAP:
885 674218 Evap.rate = Evap.monthlyEvap[0] / UCF(EVAPRATE);
886 674218 break;
887
888 90 case MONTHLY_EVAP:
889 90 Evap.rate = Evap.monthlyEvap[mon-1] / UCF(EVAPRATE);
890 90 break;
891
892 576 case TIMESERIES_EVAP:
893
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 574 times.
 576 if ( theDate >= NextEvapDate )
894 2 Evap.rate = NextEvapRate / UCF(EVAPRATE);
895 576 break;
896
897 833 case FILE_EVAP:
898 833 Evap.rate = FileValue[EVAP] / UCF(EVAPRATE);
899 833 Evap.rate *= Evap.panCoeff[mon-1];
900 833 break;
901
902 10759 case TEMPERATURE_EVAP:
903 10759 Evap.rate = FileValue[EVAP] / UCF(EVAPRATE);
904 10759 break;
905
906 ✗ default: Evap.rate = 0.0;
907 }
908
909 // --- apply climate change adjustment
910 686476 Evap.rate += Adjust.evap[mon-1];
911
912 // --- set soil recovery factor
913 686476 Evap.recoveryFactor = 1.0;
914 686476 k = Evap.recoveryPattern;
915
3/4
✓ Branch 0 taken 5760 times.
✓ Branch 1 taken 680716 times.
✓ Branch 2 taken 5760 times.
✗ Branch 3 not taken.
 686476 if ( k >= 0 && Pattern[k].type == MONTHLY_PATTERN )
916 {
917 5760 Evap.recoveryFactor = Pattern[k].factor[mon-1];
918 }
919 686476 }
920
921 //=============================================================================
922
923 686476 void setWind(DateTime theDate)
924 //
925 // Input: theDate = simulation date
926 // Output: none
927 // Purpose: sets wind speed (mph) for a specified date.
928 //
929 {
930 int yr, mon, day;
931
932
2/3
✓ Branch 0 taken 686062 times.
✓ Branch 1 taken 414 times.
✗ Branch 2 not taken.
 686476 switch ( Wind.type )
933 {
934 686062 case MONTHLY_WIND:
935 686062 datetime_decodeDate(theDate, &yr, &mon, &day);
936 686062 Wind.ws = Wind.aws[mon-1] / UCF(WINDSPEED);
937 686062 break;
938
939 414 case FILE_WIND:
940 414 Wind.ws = FileValue[WIND];
941 414 break;
942
943 ✗ default: Wind.ws = 0.0;
944 }
945 686476 }
946
947 //=============================================================================
948
949 52 void updateTempTimes(int day)
950 //
951 // Input: day = day of year
952 // Output: none
953 // Purpose: computes time of day when min/max temperatures occur.
954 // (min. temp occurs at sunrise, max. temp. at 3 hrs. < sunset)
955 //
956 {
957 double decl; // earth's declination
958 double hrang; // hour angle of sunrise/sunset
959 double arg;
960
961 52 decl = 0.40928*cos(0.017202*(172.0-day));
962 52 arg = -tan(decl)*Temp.tanAnglat;
963
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 if ( arg <= -1.0 ) arg = PI;
964
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52 times.
 52 else if ( arg >= 1.0 ) arg = 0.0;
965 52 else arg = acos(arg);
966 52 hrang = 3.8197 * arg;
967 52 Hrsr = 12.0 - hrang + Temp.dtlong;
968 52 Hrss = 12.0 + hrang + Temp.dtlong - 3.0;
969 52 Dhrdy = Hrsr - Hrss;
970 52 Dydif = 24.0 + Hrsr - Hrss;
971 52 Hrday = (Hrsr + Hrss) / 2.0;
972 52 Tave = (Tmin + Tmax) / 2.0;
973 52 Trng = (Tmax - Tmin) / 2.0;
974
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 41 times.
 52 if ( Temp.tmax == MISSING ) Trng1 = Tmax - Tmin;
975 41 else Trng1 = Temp.tmax - Tmin;
976 52 Temp.tmax = Tmax;
977 52 }
978
979 //=============================================================================
980
981 44 double getTempEvap(int day, double tave, double trng)
982 //
983 // Input: day = day of year
984 // tave = 7-day average temperature (deg F)
985 // trng = 7-day average daily temperature range (deg F)
986 // Output: returns evaporation rate in user's units (US:in/day, SI:mm/day)
987 // Purpose: uses Hargreaves method to compute daily evaporation rate
988 // from daily average temperatures and Julian day.
989 //
990 {
991 44 double a = 2.0*PI/365.0;
992 44 double ta = (tave - 32.0)*5.0/9.0; //average temperature (deg C)
993 44 double tr = trng*5.0/9.0; //temperature range (deg C)
994 44 double lamda = 2.50 - 0.002361 * ta; //latent heat of vaporization
995 44 double dr = 1.0 + 0.033*cos(a*day); //relative earth-sun distance
996 44 double phi = Temp.anglat*2.0*PI/360.0; //latitude angle (rad)
997 44 double del = 0.4093*sin(a*(284.+(double)day)); //solar declination angle (rad)
998 44 double omega = acos(-tan(phi)*tan(del)); //sunset hour angle (rad)
999 44 double ra = 37.6*dr* //extraterrestrial radiation
1000 44 (omega*sin(phi)*sin(del) +
1001 44 cos(phi)*cos(del)*sin(omega));
1002 44 double e = 0.0023*ra/lamda*sqrt(tr)*(ta+17.8); //evap. rate (mm/day)
1003
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
 44 if ( e < 0.0 ) e = 0.0;
1004
1/2
✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
 44 if ( UnitSystem == US ) e /= MMperINCH; //evap rate (in/day)
1005 44 return e;
1006 }
1007
1008 //=============================================================================
1009
1010 11 int getFileFormat()
1011 //
1012 // Input: none
1013 // Output: returns code number of climate file's format
1014 // Purpose: determines what format the climate file is in.
1015 //
1016 {
1017 11 char recdType[4] = "";
1018 11 char elemType[4] = "";
1019 11 char filler[5] = "";
1020 char staID[80];
1021 char s[80];
1022 char line[MAXLINE];
1023
1024 int y, m, d, n;
1025
1026 // --- read first line of file
1027
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 11 times.
 11 if ( fgets(line, MAXLINE, Fclimate.file) == NULL ) return UNKNOWN_FORMAT;
1028
1029 // --- check for TD3200 format
1030 11 sstrncpy(recdType, line, 3);
1031 11 sstrncpy(filler, &line[23], 4);
1032
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 10 times.
 11 if ( strcmp(recdType, "DLY") == 0 &&
1033
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 strcmp(filler, "9999") == 0 ) return TD3200;
1034
1035 // --- check for DLY0204 format
1036
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9 times.
 10 if ( strlen(line) >= 233 )
1037 {
1038 1 sstrncpy(elemType, &line[13], 3);
1039 1 n = atoi(elemType);
1040
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 1 if ( n == 1 || n == 2 || n == 151 ) return DLY0204;
1041 }
1042
1043 // --- check for USER_PREPARED format
1044 9 n = sscanf(line, "%s %d %d %d %s", staID, &y, &m, &d, s);
1045
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 2 times.
 9 if ( n == 5 ) return USER_PREPARED;
1046
1047 // --- check for GHCND format
1048
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 if ( isGhcndFormat(line) ) return GHCND;
1049
1050 ✗ return UNKNOWN_FORMAT;
1051 }
1052
1053 //=============================================================================
1054
1055 29515 void readFileLine(int *y, int *m)
1056 //
1057 // Input: none
1058 // Output: y = year
1059 // m = month
1060 // Purpose: reads year & month from next line of climate file.
1061 //
1062 {
1063 // --- read next line from climate data file
1064
2/2
✓ Branch 0 taken 29503 times.
✓ Branch 1 taken 29505 times.
 88523 while ( strlen(FileLine) == 0 )
1065 {
1066
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 29493 times.
 29503 if ( fgets(FileLine, MAXLINE, Fclimate.file) == NULL ) return;
1067
1/2
✓ Branch 0 taken 29493 times.
✗ Branch 1 not taken.
 29493 if ( FileLine[0] == '\n' ) FileLine[0] = '\0';
1068 }
1069
1070 // --- parse year & month from line
1071
4/5
✓ Branch 0 taken 29455 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 41 times.
✗ Branch 4 not taken.
 29505 switch (FileFormat)
1072 {
1073 29455 case USER_PREPARED: readUserFileLine(y, m); break;
1074 5 case TD3200: readTD3200FileLine(y,m); break;
1075 4 case DLY0204: readDLY0204FileLine(y,m); break;
1076 41 case GHCND: readGhcndFileLine(y,m); break;
1077 }
1078 }
1079
1080 //=============================================================================
1081
1082 29455 void readUserFileLine(int* y, int* m)
1083 //
1084 // Input: none
1085 // Output: y = year
1086 // m = month
1087 // Purpose: reads year & month from line of User-Prepared climate file.
1088 //
1089 {
1090 int n;
1091 char staID[80];
1092 29455 n = sscanf(FileLine, "%s %d %d", staID, y, m);
1093
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 29455 times.
 29455 if ( n < 3 )
1094 {
1095 ✗ report_writeErrorMsg(ERR_CLIMATE_FILE_READ, Fclimate.name);
1096 }
1097 29455 }
1098
1099 //=============================================================================
1100
1101 5 void readTD3200FileLine(int* y, int* m)
1102 //
1103 // Input: none
1104 // Output: y = year
1105 // m = month
1106 // Purpose: reads year & month from line of TD-3200 climate file.
1107 //
1108 {
1109 5 char recdType[4] = "";
1110 5 char year[5] = "";
1111 5 char month[3] = "";
1112
1113 // --- check for minimum number of characters
1114
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( strlen(FileLine) < 30 )
1115 {
1116 ✗ report_writeErrorMsg(ERR_CLIMATE_FILE_READ, Fclimate.name);
1117 ✗ return;
1118 }
1119
1120 // --- check for proper type of record
1121 5 sstrncpy(recdType, FileLine, 3);
1122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( strcmp(recdType, "DLY") != 0 )
1123 {
1124 ✗ report_writeErrorMsg(ERR_CLIMATE_FILE_READ, Fclimate.name);
1125 ✗ return;
1126 }
1127
1128 // --- get record's date
1129 5 sstrncpy(year, &FileLine[17], 4);
1130 5 sstrncpy(month, &FileLine[21], 2);
1131 5 *y = atoi(year);
1132 5 *m = atoi(month);
1133 }
1134
1135 //=============================================================================
1136
1137 4 void readDLY0204FileLine(int* y, int* m)
1138 //
1139 // Input: none
1140 // Output: y = year
1141 // m = month
1142 // Purpose: reads year & month from line of DLY02 or DLY04 climate file.
1143 //
1144 {
1145 4 char year[5] = "";
1146 4 char month[3] = "";
1147
1148 // --- check for minimum number of characters
1149
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( strlen(FileLine) < 16 )
1150 {
1151 ✗ report_writeErrorMsg(ERR_CLIMATE_FILE_READ, Fclimate.name);
1152 ✗ return;
1153 }
1154
1155 // --- get record's date
1156 4 sstrncpy(year, &FileLine[7], 4);
1157 4 sstrncpy(month, &FileLine[11], 2);
1158 4 *y = atoi(year);
1159 4 *m = atoi(month);
1160 }
1161
1162 //=============================================================================
1163
1164 12 void readFileValues()
1165 //
1166 // Input: none
1167 // Output: none
1168 // Purpose: reads next month's worth of data from climate file.
1169 //
1170 {
1171 int i, j;
1172 int y, m;
1173
1174 // --- initialize FileData array to missing values
1175
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 12 times.
 60 for ( i=0; i<MAXCLIMATEVARS; i++)
1176 {
1177
2/2
✓ Branch 0 taken 1536 times.
✓ Branch 1 taken 48 times.
 1584 for (j=0; j<MAXDAYSPERMONTH; j++) FileData[i][j] = MISSING;
1178 }
1179
1180
1/2
✓ Branch 0 taken 158 times.
✗ Branch 1 not taken.
 158 while ( !ErrorCode )
1181 {
1182 // --- return when date on line is after current file date
1183
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 148 times.
 160 if ( feof(Fclimate.file) ) return;
1184 148 readFileLine(&y, &m);
1185
4/4
✓ Branch 0 taken 147 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 146 times.
 148 if ( y > FileYear || m > FileMonth ) return;
1186
1187 // --- parse climate values from file line
1188
4/5
✓ Branch 0 taken 98 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 39 times.
✗ Branch 4 not taken.
 146 switch (FileFormat)
1189 {
1190 98 case USER_PREPARED: parseUserFileLine(); break;
1191 5 case TD3200: parseTD3200FileLine(); break;
1192 4 case DLY0204: parseDLY0204FileLine(); break;
1193 39 case GHCND: parseGhcndFileLine(); break;
1194 }
1195 146 sstrncpy(FileLine, "", 0);
1196 }
1197 }
1198
1199 //=============================================================================
1200
1201 98 void parseUserFileLine()
1202 //
1203 // Input: none
1204 // Output: none
1205 // Purpose: parses climate variable values from a line of a user-prepared
1206 // climate file.
1207 //
1208 {
1209 int n;
1210 int y, m, d;
1211 char staID[80];
1212 98 char s0[80] = "";
1213 98 char s1[80] = "";
1214 98 char s2[80] = "";
1215 98 char s3[80] = "";
1216 double x;
1217
1218 // --- read day, Tmax, Tmin, Evap, & Wind from file line
1219 98 n = sscanf(FileLine, "%s %d %d %d %s %s %s %s",
1220 staID, &y, &m, &d, s0, s1, s2, s3);
1221
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 92 times.
 98 if ( n < 4 ) return;
1222
2/4
✓ Branch 0 taken 92 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 92 times.
 92 if ( d < 1 || d > 31 ) return;
1223
1224 // --- process TMAX
1225
2/4
✓ Branch 0 taken 92 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 92 times.
✗ Branch 3 not taken.
 92 if ( strlen(s0) > 0 && *s0 != '*' )
1226 {
1227 92 x = atof(s0);
1228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
 92 if ( UnitSystem == SI ) x = 9./5.*x + 32.0;
1229 92 FileData[TMAX][d] = x;
1230 }
1231
1232 // --- process TMIN
1233
2/4
✓ Branch 0 taken 92 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 92 times.
✗ Branch 3 not taken.
 92 if ( strlen(s1) > 0 && *s1 != '*' )
1234 {
1235 92 x = atof(s1);
1236
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
 92 if ( UnitSystem == SI ) x = 9./5.*x + 32.0;
1237 92 FileData[TMIN][d] = x;
1238 }
1239
1240 // --- process EVAP
1241
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 92 if ( strlen(s2) > 0 && *s2 != '*' ) FileData[EVAP][d] = atof(s2);
1242
1243 // --- process WIND
1244
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 92 if ( strlen(s3) > 0 && *s3 != '*' ) FileData[WIND][d] = atof(s3);
1245 }
1246
1247 //=============================================================================
1248
1249 5 void parseTD3200FileLine()
1250 //
1251 // Input: none
1252 // Output: none
1253 // Purpose: parses climate variable values from a line of a TD3200 file.
1254 //
1255 {
1256 int i;
1257 5 char param[5] = "";
1258
1259 // --- parse parameter name
1260 5 sstrncpy(param, &FileLine[11], 4);
1261
1262 // --- see if parameter is temperature, evaporation or wind speed
1263
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
 25 for (i=0; i<MAXCLIMATEVARS; i++)
1264 {
1265
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 17 times.
 20 if (strcmp(param, ClimateVarWords[i]) == 0 ) setTD3200FileValues(i);
1266 }
1267 5 }
1268
1269 //=============================================================================
1270
1271 3 void setTD3200FileValues(int i)
1272 //
1273 // Input: i = climate variable code
1274 // Output: none
1275 // Purpose: reads month worth of values for climate variable from TD-3200 file.
1276 //
1277 {
1278 3 char valCount[4] = "";
1279 3 char day[3] = "";
1280 3 char sign[2] = "";
1281 3 char value[6] = "";
1282 3 char flag2[2] = "";
1283 double x;
1284 int nValues;
1285 int j, k, d;
1286
1287 // --- parse number of days with data from cols. 27-29 of file line
1288 3 sstrncpy(valCount, &FileLine[27], 3);
1289 3 nValues = atoi(valCount);
1290
1291 // --- check for enough characters on line
1292
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if ( (int)strlen(FileLine) >= 12*nValues + 30 )
1293 {
1294 // --- for each day's value
1295
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 3 times.
 18 for (j=0; j<nValues; j++)
1296 {
1297 // --- parse day, value & flag from file line
1298 15 k = 30 + j*12;
1299 15 sstrncpy(day, &FileLine[k], 2);
1300 15 sstrncpy(sign, &FileLine[k+4], 1);
1301 15 sstrncpy(value, &FileLine[k+5], 5);
1302 15 sstrncpy(flag2, &FileLine[k+11], 1);
1303
1304 // --- if value is valid then store it in FileData array
1305 15 d = atoi(day);
1306
1/2
✓ Branch 0 taken 15 times.
✗ Branch 1 not taken.
 15 if ( strcmp(value, "99999") != 0
1307
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 15 && ( flag2[0] == '0' || flag2[0] == '1')
1308
1/2
✓ Branch 0 taken 15 times.
✗ Branch 1 not taken.
 15 && d > 0
1309
1/2
✓ Branch 0 taken 15 times.
✗ Branch 1 not taken.
 15 && d <= 31 )
1310 {
1311 // --- convert from string value to numerical value
1312 15 x = atof(value);
1313
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 if ( sign[0] == '-' ) x = -x;
1314
1315 // --- convert evaporation from hundreths of inches
1316
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10 times.
 15 if ( i == EVAP )
1317 {
1318 5 x /= 100.0;
1319
1320 // --- convert to mm if using SI units
1321
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( UnitSystem == SI ) x *= MMperINCH;
1322 }
1323
1324 // --- convert wind speed from miles/day to miles/hour
1325
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 if ( i == WIND ) x /= 24.0;
1326
1327 // --- store value
1328 15 FileData[i][d] = x;
1329 }
1330 }
1331 }
1332 3 }
1333
1334 //=============================================================================
1335
1336 4 void parseDLY0204FileLine()
1337 //
1338 // Input: none
1339 // Output: none
1340 // Purpose: parses a month's worth of climate variable values from a line of
1341 // a DLY02 or DLY04 climate file.
1342 //
1343 {
1344 int j, k, p;
1345 4 char param[4] = "";
1346 4 char sign[2] = "";
1347 4 char value[6] = "";
1348 4 char code[2] = "";
1349 double x;
1350
1351 // --- parse parameter name
1352 4 sstrncpy(param, &FileLine[13], 3);
1353
1354 // --- see if parameter is min or max temperature
1355 4 p = atoi(param);
1356
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ( p == 1 ) p = TMAX;
1357
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
 3 else if ( p == 2 ) p = TMIN;
1358
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 else if ( p == 151 ) p = EVAP;
1359 1 else return;
1360
1361 // --- check for 233 characters on line
1362
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ( strlen(FileLine) < 233 ) return;
1363
1364 // --- for each of 31 days
1365 3 k = 16;
1366
2/2
✓ Branch 0 taken 93 times.
✓ Branch 1 taken 3 times.
 96 for (j=1; j<=31; j++)
1367 {
1368 // --- parse value & flag from file line
1369 93 sstrncpy(sign, &FileLine[k], 1);
1370 93 sstrncpy(value, &FileLine[k+1], 5);
1371 93 sstrncpy(code, &FileLine[k+6], 1);
1372 93 k += 7;
1373
1374 // --- if value is valid then store it in FileData array
1375
1376
2/4
✓ Branch 0 taken 93 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 93 times.
✗ Branch 3 not taken.
 93 if ( strcmp(value, "99999") != 0 && strcmp(value, " ") != 0 )
1377 {
1378
2/3
✓ Branch 0 taken 62 times.
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 93 switch (p)
1379 {
1380 62 case TMAX:
1381 case TMIN:
1382 // --- convert from integer tenths of a degree C to degrees F
1383 62 x = atof(value) / 10.0;
1384
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 62 times.
 62 if ( sign[0] == '-' ) x = -x;
1385 62 x = 9./5.*x + 32.0;
1386 62 break;
1387 31 case EVAP:
1388 // --- convert from 0.1 mm to inches or mm
1389 31 x = atof(value) / 10.0;
1390
1/2
✓ Branch 0 taken 31 times.
✗ Branch 1 not taken.
 31 if ( UnitSystem == US ) x /= MMperINCH;
1391 31 break;
1392 ✗ default: return;
1393 }
1394 93 FileData[p][j] = x;
1395 }
1396 }
1397 }
1398
1399 //=============================================================================
1400
1401 2 int isGhcndFormat(char* line)
1402 //
1403 // Input: line = first line of text from a climate file
1404 // Output: returns TRUE if climate file is in NCDC GHCN Daily format.
1405 // Purpose: Checks if a climate file is in the NCDC GHCN Daily format
1406 // and determines the position of each climate variable field.
1407 //
1408 {
1409 int i;
1410 char* ptr;
1411
1412 // --- find starting position of the DATE field
1413 2 ptr = strstr(line, "DATE");
1414
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ptr == NULL ) return FALSE;
1415 2 FileDateFieldPos = (int)(ptr - line);
1416
1417 // --- initialize starting position of each data field
1418
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for ( i = TMIN; i <= WIND; i++) FileFieldPos[i] = -1;
1419
1420 // --- find starting position of each climate variable's data field
1421 2 ptr = strstr(line, "TMIN");
1422
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( ptr ) FileFieldPos[TMIN] = (int)(ptr - line);
1423 2 ptr = strstr(line, "TMAX");
1424
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( ptr ) FileFieldPos[TMAX] = (int)(ptr - line);
1425 2 ptr = strstr(line, "EVAP");
1426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ptr ) FileFieldPos[EVAP] = (int)(ptr - line);
1427
1428 // --- WIND can either be daily movement or average speed
1429 2 FileWindType = WDMV;
1430 2 ptr = strstr(line, "WDMV");
1431
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( ptr == NULL )
1432 {
1433 2 FileWindType = AWND;
1434 2 ptr = strstr(line, "AWND");
1435 }
1436
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ptr ) FileFieldPos[WIND] = (int)(ptr - line);
1437
1438 // --- check if at least one climate variable was found
1439
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 for (i = TMIN; i <= WIND; i++) if (FileFieldPos[i] >= 0 ) return TRUE;
1440 ✗ return FALSE;
1441 }
1442
1443 //=============================================================================
1444
1445 41 void readGhcndFileLine(int* y, int* m)
1446 //
1447 // Input: none
1448 // Output: y = year
1449 // m = month
1450 // Purpose: reads year & month from line of a NCDC GHCN Daily climate file.
1451 //
1452 {
1453 41 int n = sscanf(&FileLine[FileDateFieldPos], "%4d%2d", y, m);
1454
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 39 times.
 41 if ( n != 2 )
1455 {
1456 2 *y = -99999;
1457 2 *m = -99999;
1458 }
1459 41 }
1460
1461 //=============================================================================
1462
1463 39 void parseGhcndFileLine()
1464 //
1465 // Input: none
1466 // Output: none
1467 // Purpose: parses a line of a NCDC GHCN Daily file for daily
1468 // values of max/min temperature, pan evaporation and
1469 // wind speed.
1470 //
1471 {
1472 int y, m, d, n, i;
1473 double v;
1474
1475 // --- parse day of month from date field
1476 39 n = sscanf(&FileLine[FileDateFieldPos], "%4d%2d%2d", &y, &m, &d);
1477
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 37 times.
 39 if ( n < 3 ) return;
1478
2/4
✓ Branch 0 taken 37 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 37 times.
 37 if ( d < 1 || d > 31 ) return;
1479
1480 // --- parse climate variables
1481
2/2
✓ Branch 0 taken 148 times.
✓ Branch 1 taken 37 times.
 185 for (i = TMIN; i <= WIND; i++)
1482 {
1483
2/2
✓ Branch 0 taken 74 times.
✓ Branch 1 taken 74 times.
 148 if ( FileFieldPos[i] >= 0 )
1484 {
1485
1/2
✓ Branch 0 taken 74 times.
✗ Branch 1 not taken.
 74 if ( sscanf(&FileLine[FileFieldPos[i]], "%8lf", &v) > 0 )
1486 {
1487
1/2
✓ Branch 0 taken 74 times.
✗ Branch 1 not taken.
 74 if ( fabs(v) < 9999. )
1488 74 FileData[i][d] = convertGhcndValue(i, v);
1489 }
1490 }
1491 }
1492 }
1493
1494 //=============================================================================
1495
1496 74 double convertGhcndValue(int var, double v)
1497 //
1498 // Input: var = climate variable code
1499 // v = climate variable value
1500 // Output: climate variable value in SWMM's internal units
1501 // Purpose: converts a climate variable value read from a NCDC GHCN Daily file
1502 // to SWMM's internal units.
1503 //
1504 {
1505
1/4
✓ Branch 0 taken 74 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 74 switch (var)
1506 {
1507 74 case TMIN:
1508 case TMAX:
1509
1/3
✓ Branch 0 taken 74 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 74 switch (FileTempUnits)
1510 {
1511 74 case DEG_C10: // tenths deg. C ==> deg. F
1512 74 return v / 10. * 9.0 / 5.0 + 32.0;
1513
1514 ✗ case DEG_C: // deg. C ==> deg. F
1515 ✗ return v * 9.0 / 5.0 + 32.0;
1516
1517 ✗ default: // deg. F
1518 ✗ return v;
1519 }
1520 ✗ case EVAP:
1521 ✗ switch (FileTempUnits)
1522 {
1523 ✗ case DEG_C10: // tenths mm ==> inches or mm
1524 ✗ v /= 10.;
1525 ✗ if (UnitSystem == US) v /= MMperINCH;
1526 ✗ return v;
1527
1528 ✗ case DEG_C: // mm ==> inches or mm
1529 ✗ if (UnitSystem == US) v /= MMperINCH;
1530 ✗ return v;
1531
1532 ✗ default: // inches ==> inches or mm
1533 ✗ if (UnitSystem == SI) v *= MMperINCH;
1534 ✗ return v;
1535 }
1536 ✗ case WIND:
1537 ✗ switch (FileTempUnits)
1538 {
1539 ✗ case DEG_C10:
1540 // km/day ==> miles/hr
1541 ✗ if (FileWindType == WDMV)
1542 ✗ return v * 0.62137 / 24.;
1543 // tenths m/s ==> miles/hr
1544 else
1545 ✗ return v / 10. / 1000. * 0.62137 * 3600.;
1546 ✗ case DEG_C:
1547 // km/day ==> miles/hr
1548 ✗ if (FileWindType == WDMV)
1549 ✗ return v * 0.62137 / 24.;
1550 // m/s ==> miles/hr
1551 else
1552 ✗ return v / 1000. * 0.62137 * 3600.;
1553
1554 ✗ default:
1555 // miles ==> miles/hr
1556 ✗ if (FileWindType == WDMV)
1557 ✗ return v / 24.;
1558 // miles/hr
1559 else
1560 ✗ return v;
1561 }
1562 ✗ default:
1563 ✗ return v;
1564 }
1565 }
1566
1567 //=============================================================================
1568
1569 44 void updateTempMoveAve(double tmin, double tmax)
1570 //
1571 // Input: tmin = minimum daily temperature (deg F)
1572 // tmax = maximum daily temperature (deg F)
1573 // Output: none
1574 // Purpose: updates moving averages of average daily temperature
1575 // and daily temperature range stored in structure Tma.
1576 //
1577 {
1578 double ta, // new day's average temperature (deg F)
1579 tr; // new day's temperature range (deg F)
1580 44 int kount = Tma.count;
1581 44 double count = kount;
1582
1583 // --- find ta and tr from new day's min and max temperature
1584 44 ta = (tmin + tmax) / 2.0;
1585 44 tr = fabs(tmax - tmin);
1586
1587 // --- if the array used to store previous days' temperatures is full
1588
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 19 times.
 44 if ( kount == Tma.maxCount )
1589 {
1590 // --- update the moving averages with the new day's value
1591 25 Tma.tAve = (Tma.tAve * count + ta - Tma.ta[Tma.front]) / count;
1592 25 Tma.tRng = (Tma.tRng * count + tr - Tma.tr[Tma.front]) / count;
1593
1594 // --- replace the values at the front of the moving average window
1595 25 Tma.ta[Tma.front] = ta;
1596 25 Tma.tr[Tma.front] = tr;
1597
1598 // --- move the front one position forward
1599 25 Tma.front++;
1600
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 22 times.
 25 if ( Tma.front == count ) Tma.front = 0;
1601 }
1602
1603 // --- array of previous day's values not full (at start of simulation)
1604 else
1605 {
1606 // --- find new moving averages by adding new values to previous ones
1607 19 Tma.tAve = (Tma.tAve * count + ta) / (count + 1);
1608 19 Tma.tRng = (Tma.tRng * count + tr) / (count + 1);
1609
1610 // --- save new day's values
1611 19 Tma.ta[Tma.front] = ta;
1612 19 Tma.tr[Tma.front] = tr;
1613
1614 // --- increment count and front of moving average window
1615 19 Tma.count++;
1616 19 Tma.front++;
1617
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 18 times.
 19 if ( Tma.count == Tma.maxCount ) Tma.front = 0;
1618 }
1619 44 }
1620