1 #include <stdlib.h>
2 #include <string.h>
3 #include <time.h>
4 #include <grass/raster.h>
5 #include <grass/glocale.h>
6 #include "global.h"
7
get_line(FILE * fp,char * buffer)8 void get_line(FILE * fp, char *buffer)
9 {
10 char *str;
11
12 buffer[0] = 0;
13 fscanf(fp, "%[^\n]", buffer);
14 getc(fp);
15
16 if ((str = (char *)strchr(buffer, '#')))
17 *str = 0;
18 }
19
read_input(void)20 void read_input(void)
21 {
22 char buf[1024];
23 FILE *fp;
24 int i;
25
26 /* Read parameters file */
27 if ((fp = fopen(file.params, "r")) == NULL)
28 G_fatal_error(_("Unable to open input file <%s>"), file.params);
29
30 for (; !feof(fp);) {
31 get_line(fp, buf);
32 i = strlen(buf) - 1;
33 for (; i >= 0; i--) {
34 if (buf[i] != ' ' && buf[i] != '\t') {
35 buf[i + 1] = 0;
36 break;
37 }
38 }
39 if (i >= 0)
40 break;
41 }
42 params.name = G_store(buf);
43
44 for (; !feof(fp);) {
45 get_line(fp, buf);
46 if (sscanf(buf, "%lf", &(params.A)) == 1)
47 break;
48 }
49 for (; !feof(fp);) {
50 get_line(fp, buf);
51 if (sscanf(buf, "%lf", &(params.qs0)) == 1)
52 break;
53 }
54 if (params.qs0 == 0.0) {
55 fclose(fp);
56 G_fatal_error(_("%s cannot be 0"), "parameters.qs0");
57 }
58 for (; !feof(fp);) {
59 get_line(fp, buf);
60 if (sscanf(buf, "%lf", &(params.lnTe)) == 1)
61 break;
62 }
63 for (; !feof(fp);) {
64 get_line(fp, buf);
65 if (sscanf(buf, "%lf", &(params.m)) == 1)
66 break;
67 }
68 for (; !feof(fp);) {
69 get_line(fp, buf);
70 if (sscanf(buf, "%lf", &(params.Sr0)) == 1)
71 break;
72 }
73 for (; !feof(fp);) {
74 get_line(fp, buf);
75 if (sscanf(buf, "%lf", &(params.Srmax)) == 1)
76 break;
77 }
78 for (; !feof(fp);) {
79 get_line(fp, buf);
80 if (sscanf(buf, "%lf", &(params.td)) == 1)
81 break;
82 }
83 for (; !feof(fp);) {
84 get_line(fp, buf);
85 if (sscanf(buf, "%lf", &(params.vch)) == 1)
86 break;
87 }
88 for (; !feof(fp);) {
89 get_line(fp, buf);
90 if (sscanf(buf, "%lf", &(params.vr)) == 1)
91 break;
92 }
93 for (; !feof(fp);) {
94 get_line(fp, buf);
95 if (sscanf(buf, "%d", &(params.infex)) == 1)
96 break;
97 }
98 for (; !feof(fp);) {
99 get_line(fp, buf);
100 if (sscanf(buf, "%lf", &(params.K0)) == 1)
101 break;
102 }
103 for (; !feof(fp);) {
104 get_line(fp, buf);
105 if (sscanf(buf, "%lf", &(params.psi)) == 1)
106 break;
107 }
108 for (; !feof(fp);) {
109 get_line(fp, buf);
110 if (sscanf(buf, "%lf", &(params.dtheta)) == 1)
111 break;
112 }
113
114 params.d = NULL;
115 params.Ad = NULL;
116
117 for (i = 0; !feof(fp);) {
118 double d;
119 double Ad_r;
120
121 get_line(fp, buf);
122 if (sscanf(buf, "%lf %lf", &d, &Ad_r) == 2) {
123 params.d = (double *)G_realloc(params.d, (i + 1) * sizeof(double));
124 params.Ad = (double *)G_realloc(params.Ad,
125 (i + 1) * sizeof(double));
126 params.d[i] = d;
127 params.Ad[i++] = Ad_r * params.A;
128 }
129 }
130
131 params.nch = i;
132 fclose(fp);
133
134 /* Read topographic index statistics file */
135 if ((fp = fopen(file.topidxstats, "r")) == NULL)
136 G_fatal_error(_("Unable to open input file <%s>"), file.topidxstats);
137
138 topidxstats.atb = NULL;
139 topidxstats.Aatb_r = NULL;
140
141 for (i = 0; !feof(fp);) {
142 double atb;
143 double Aatb_r;
144
145 get_line(fp, buf);
146 if (sscanf(buf, "%lf %lf", &atb, &Aatb_r) == 2) {
147 topidxstats.atb = (double *)G_realloc(topidxstats.atb,
148 (i + 1) * sizeof(double));
149 topidxstats.Aatb_r = (double *)G_realloc(topidxstats.Aatb_r,
150 (i + 1) * sizeof(double));
151 topidxstats.atb[i] = atb;
152 topidxstats.Aatb_r[i++] = Aatb_r;
153 }
154 }
155
156 misc.ntopidxclasses = i;
157 fclose(fp);
158
159 /* Read input file */
160 if ((fp = fopen(file.input, "r")) == NULL)
161 G_fatal_error(_("Unable to open input file <%s>"), file.input);
162
163 for (; !feof(fp);) {
164 get_line(fp, buf);
165 if (sscanf(buf, "%lf", &(input.dt)) == 1)
166 break;
167 }
168
169 input.R = NULL;
170 input.Ep = NULL;
171
172 for (i = 0; !feof(fp);) {
173 double R;
174 double Ep;
175
176 get_line(fp, buf);
177 if (sscanf(buf, "%lf %lf", &R, &Ep) == 2) {
178 input.R = (double *)G_realloc(input.R, (i + 1) * sizeof(double));
179 input.Ep = (double *)G_realloc(input.Ep, (i + 1) * sizeof(double));
180 input.R[i] = R;
181 input.Ep[i++] = Ep;
182 }
183 }
184
185 input.ntimesteps = i;
186 fclose(fp);
187
188 if (!(misc.timestep > 0 && misc.timestep < input.ntimesteps + 1))
189 misc.timestep = 0;
190 if (!(misc.topidxclass > 0 && misc.topidxclass < misc.ntopidxclasses + 1))
191 misc.topidxclass = 0;
192 }
193
write_output(void)194 void write_output(void)
195 {
196 FILE *fp;
197 time_t tloc;
198 struct tm *ltime;
199 int st, et, si, ei;
200 int i, j;
201
202 time(&tloc);
203 ltime = localtime(&tloc);
204
205 ltime->tm_year += 1900;
206 ltime->tm_mon++;
207
208 if ((fp = fopen(file.output, "w")) == NULL)
209 G_fatal_error(_("Unable to create output file <%s>"), file.output);
210
211 fprintf(fp, "# r.topmodel output file for %s\n", params.name);
212 fprintf(fp, "# Run time: %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
213 ltime->tm_year, ltime->tm_mon, ltime->tm_mday,
214 ltime->tm_hour, ltime->tm_min, ltime->tm_sec);
215 fprintf(fp, "#\n");
216 fprintf(fp, "# Qt_peak [m^3/timestep]: Peak total flow\n");
217 fprintf(fp, "# tt_peak [timestep]: Peak time for total flow\n");
218 fprintf(fp, "# Qt_mean [m^3/timestep]: Mean total flow\n");
219 fprintf(fp, "# lnTe [ln(m^2/timestep)]: ln of the average transmissivity at the soil surface\n");
220 fprintf(fp, "# vch [m/timestep]: Main channel routing velocity\n");
221 fprintf(fp, "# vr [m/timestep]: Internal subcatchment routing velocity\n");
222 fprintf(fp, "# lambda [ln(m^2)]: Average topographic index\n");
223 fprintf(fp, "# qss [m/timestep]: Saturated subsurface flow per unit area\n");
224 fprintf(fp, "# qs0 [m/timestep]: Initial subsurface flow per unit area\n");
225 fprintf(fp, "# ntopidxclasses: Number of topographic index classes\n");
226 fprintf(fp, "# dt [h]: Time step\n");
227 fprintf(fp, "# ntimesteps: Number of time steps\n");
228 fprintf(fp, "# nch: Number of channel segments\n");
229 fprintf(fp, "# delay [timestep]: Routing delay in the main channel\n");
230 fprintf(fp, "# tc [timestep]: Time of concentration\n");
231 fprintf(fp, "# tcsub [timestep]: Time of concentration in the subcatchment\n");
232 fprintf(fp, "#\n");
233 fprintf(fp, "# tch [timestep]: Routing time to the catchment outlet\n");
234 fprintf(fp, "# Ad [m^2]: Difference in the contribution area\n");
235 fprintf(fp, "# Qt [m^3/timestep]: Total flow\n");
236 fprintf(fp, "# qt [m/timestep]: Total flow per unit area\n");
237 fprintf(fp, "# qo [m/timestep]: Saturated overland flow per unit area\n");
238 fprintf(fp, "# qs [m/timestep]: Subsurface flow per unit area\n");
239 fprintf(fp, "# qv [m/timestep]: Vertical drainage flux from unsaturated zone\n");
240 fprintf(fp, "# S_mean [m]: Mean saturation deficit\n");
241 if (params.infex) {
242 fprintf(fp, "# f [m/timestep]: Infiltration rate\n");
243 fprintf(fp, "# fex [m/timestep]: Infiltration excess runoff\n");
244 }
245
246 if (misc.timestep || misc.topidxclass) {
247 fprintf(fp, "#\n");
248 fprintf(fp, "# Srz [m]: Root zone storage deficit\n");
249 fprintf(fp, "# Suz [m]: Unsaturated zone storage (gravity drainage)\n");
250 fprintf(fp, "# S [m]: Local saturated zone deficit due to gravity drainage\n");
251 fprintf(fp, "# Ea [m/timestep]: Actual evapotranspiration\n");
252 fprintf(fp, "# ex [m/timestep]: Excess flow from fully saturated area per unit area\n");
253 }
254 fprintf(fp, "\n");
255
256 fprintf(fp, "Qt_peak: %10.3e\n", misc.Qt_peak);
257 fprintf(fp, "tt_peak: %10d\n", misc.tt_peak);
258 fprintf(fp, "Qt_mean: %10.3e\n", misc.Qt_mean);
259 fprintf(fp, "lnTe: %10.3e\n", misc.lnTe);
260 fprintf(fp, "vch: %10.3e\n", misc.vch);
261 fprintf(fp, "vr: %10.3e\n", misc.vr);
262 fprintf(fp, "lambda: %10.3e\n", misc.lambda);
263 fprintf(fp, "qss: %10.3e\n", misc.qss);
264 fprintf(fp, "qs0: %10.3e\n", misc.qs0);
265 fprintf(fp, "ntopidxclasses: %10d\n", misc.ntopidxclasses);
266 fprintf(fp, "dt: %10.3e\n", input.dt);
267 fprintf(fp, "ntimesteps: %10d\n", input.ntimesteps);
268 fprintf(fp, "nch: %10d\n", params.nch);
269 fprintf(fp, "delay: %10d\n", misc.delay);
270 fprintf(fp, "tc: %10d\n", misc.tc);
271 fprintf(fp, "tcsub: %10d\n", misc.tcsub);
272 fprintf(fp, "\n");
273
274 fprintf(fp, "%10s\n", "tch");
275 for (i = 0; i < params.nch; i++)
276 fprintf(fp, "%10.3e\n", misc.tch[i]);
277 fprintf(fp, "\n");
278
279 fprintf(fp, "%10s\n", "Ad");
280 for (i = 0; i < misc.tcsub; i++)
281 fprintf(fp, "%10.3e\n", misc.Ad[i]);
282 fprintf(fp, "\n");
283
284 st = et = si = ei = 0;
285 if (misc.timestep || misc.topidxclass) {
286 if (misc.timestep) {
287 st = misc.timestep - 1;
288 et = misc.timestep;
289 }
290 else {
291 st = 0;
292 et = input.ntimesteps;
293 }
294
295 if (misc.topidxclass) {
296 si = misc.topidxclass - 1;
297 ei = misc.topidxclass;
298 }
299 else {
300 si = 0;
301 ei = misc.ntopidxclasses;
302 }
303 }
304
305 fprintf(fp, "%10s %10s %10s %10s %10s %10s %10s",
306 "timestep", "Qt", "qt", "qo", "qs", "qv", "S_mean");
307 if (params.infex)
308 fprintf(fp, " %10s %10s", "f", "fex");
309 fprintf(fp, "\n");
310
311 for (i = 0; i < input.ntimesteps; i++) {
312 fprintf(fp, "%10d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e",
313 i + 1, misc.Qt[i],
314 misc.qt[i][misc.ntopidxclasses],
315 misc.qo[i][misc.ntopidxclasses], misc.qs[i],
316 misc.qv[i][misc.ntopidxclasses], misc.S_mean[i]);
317 if (params.infex)
318 fprintf(fp, " %10.3e %10.3e", misc.f[i], misc.fex[i]);
319 fprintf(fp, "\n");
320 }
321
322 if (misc.timestep || misc.topidxclass) {
323 fprintf(fp, "\n");
324 fprintf(fp, "For ");
325 if (misc.timestep)
326 fprintf(fp, "timestep: %5d", misc.timestep);
327 if (misc.timestep && misc.topidxclass)
328 fprintf(fp, ", ");
329 if (misc.topidxclass)
330 fprintf(fp, "topidxclass: %5d", misc.topidxclass);
331 fprintf(fp, "\n");
332
333 if (misc.timestep && !misc.topidxclass) {
334 fprintf(fp, "%10s ", "topidxclass");
335 }
336 else if (misc.topidxclass && !misc.timestep) {
337 fprintf(fp, "%10s ", "timestep");
338 }
339
340 fprintf(fp, "%10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
341 "qt", "qo", "qs", "qv", "Srz", "Suz", "S", "Ea", "ex");
342
343 for (i = st; i < et; i++)
344 for (j = si; j < ei; j++) {
345 if (misc.timestep && !misc.topidxclass) {
346 fprintf(fp, "%10d ", j + 1);
347 }
348 else if (misc.topidxclass && !misc.timestep) {
349 fprintf(fp, "%10d ", i + 1);
350 }
351
352 fprintf(fp, "%10.3e %10.3e %10.3e "
353 "%10.3e %10.3e %10.3e "
354 "%10.3e %10.3e %10.3e\n",
355 misc.qt[i][j], misc.qo[i][j],
356 misc.qs[i], misc.qv[i][j],
357 misc.Srz[i][j], misc.Suz[i][j],
358 misc.S[i][j], misc.Ea[i][j], misc.ex[i][j]);
359 }
360 }
361
362 fclose(fp);
363 }
364