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