Used to do too much preprocessing. Now, it optionally takes the masked topidx

map only for generating the topographic index statistics file. This change
reduces user confusion and simplifies GUI.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@58073 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.topmodel/file_io.c

@@ -1,8 +1,11 @@
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
+#include <grass/glocale.h>
 #include "global.h"
 
-
 void get_line(FILE * fp, char *buffer)
 {
     char *str;
@@ -15,12 +18,10 @@ void get_line(FILE * fp, char *buffer)
     if ((str = (char *)strchr(buffer, '#')))
 	*str = 0;
 
-
     return;
 }
 
-
-void read_inputs(void)
+void read_input(void)
 {
     char buf[1024];
     FILE *fp;
@@ -29,42 +30,54 @@ void read_inputs(void)
 
 
     /* Read topographic index statistics file */
-    fp = fopen(file.idxstats, "r");
-    idxstats.atb = (double *)G_malloc(misc.nidxclass * sizeof(double));
-    idxstats.Aatb_r = (double *)G_malloc(misc.nidxclass * sizeof(double));
+    if ((fp = fopen(file.topidxstats, "r")) == NULL)
+	G_fatal_error(_("Unable to open input file <%s>"), file.topidxstats);
+
+    topidxstats.atb = NULL;
+    topidxstats.Aatb_r = NULL;
+    misc.ncells = 0;
 
-    misc.ncell = 0;
+    for (i = 0, j = 0; !feof(fp);) {
+	double atb;
+	double Aatb_r;
 
-    for (i = 0; i < misc.nidxclass && !feof(fp);) {
 	get_line(fp, buf);
 
-	if (sscanf(buf, "%lf %lf",
-		   &(idxstats.atb[i]), &(idxstats.Aatb_r[i])) == 2)
-	    misc.ncell += (int)idxstats.Aatb_r[i++];
+	if (sscanf(buf, "%lf %lf", &atb, &Aatb_r) == 2) {
+	    topidxstats.atb = (double *)G_realloc(topidxstats.atb,
+			    (j + 1) * sizeof(double));
+	    topidxstats.Aatb_r = (double *)G_realloc(topidxstats.Aatb_r,
+			    (j + 1) * sizeof(double));
+	    topidxstats.atb[j] = atb;
+	    topidxstats.Aatb_r[j] = Aatb_r;
+	    misc.ncells += (int)topidxstats.Aatb_r[j++];
+	}
     }
 
-    misc.nidxclass = i;
+    misc.ntopidxclasses = j;
+
     fclose(fp);
 
-    for (i = 0; i < misc.nidxclass; i++)
-	idxstats.Aatb_r[i] /= (double)misc.ncell;
-
-    for (i = 0; i < misc.nidxclass; i++) {
-	for (j = i; j < misc.nidxclass; j++) {
-	    if (idxstats.atb[i] < idxstats.atb[j]) {
-		x = idxstats.atb[i];
-		idxstats.atb[i] = idxstats.atb[j];
-		idxstats.atb[j] = x;
-		x = idxstats.Aatb_r[i];
-		idxstats.Aatb_r[i] = idxstats.Aatb_r[j];
-		idxstats.Aatb_r[j] = x;
+    for (i = 0; i < misc.ntopidxclasses; i++)
+	topidxstats.Aatb_r[i] /= (double)misc.ncells;
+
+    for (i = 0; i < misc.ntopidxclasses; i++) {
+	for (j = i; j < misc.ntopidxclasses; j++) {
+	    if (topidxstats.atb[i] < topidxstats.atb[j]) {
+		x = topidxstats.atb[i];
+		topidxstats.atb[i] = topidxstats.atb[j];
+		topidxstats.atb[j] = x;
+		x = topidxstats.Aatb_r[i];
+		topidxstats.Aatb_r[i] = topidxstats.Aatb_r[j];
+		topidxstats.Aatb_r[j] = x;
 	    }
 	}
     }
 
 
     /* Read parameters file */
-    fp = fopen(file.params, "r");
+    if ((fp = fopen(file.params, "r")) == NULL)
+	G_fatal_error(_("Unable to open input file <%s>"), file.params);
 
     for (; !feof(fp);) {
 	get_line(fp, buf);
@@ -136,58 +149,57 @@ void read_inputs(void)
 
 
     /* Read input file */
-    fp = fopen(file.input, "r");
+    if ((fp = fopen(file.input, "r")) == NULL)
+	G_fatal_error(_("Unable to open input file <%s>"), file.input);
 
     for (; !feof(fp);) {
 	get_line(fp, buf);
 
-	if (sscanf(buf, "%d %lf", &(input.ntimestep), &(input.dt)) == 2)
+	if (sscanf(buf, "%d %lf", &(input.ntimesteps), &(input.dt)) == 2)
 	    break;
     }
 
-    input.R = (double *)G_malloc(input.ntimestep * sizeof(double));
-    input.Ep = (double *)G_malloc(input.ntimestep * sizeof(double));
+    input.R = (double *)G_malloc(input.ntimesteps * sizeof(double));
+    input.Ep = (double *)G_malloc(input.ntimesteps * sizeof(double));
 
-    for (i = 0; i < input.ntimestep && !feof(fp);) {
+    for (i = 0; i < input.ntimesteps && !feof(fp);) {
 	get_line(fp, buf);
 
 	if (sscanf(buf, "%lf %lf", &(input.R[i]), &(input.Ep[i])) == 2)
 	    i++;
     }
 
-    input.ntimestep = i;
+    input.ntimesteps = i;
     fclose(fp);
 
-
     /* Read Qobs file */
-    if (file.Qobs) {
-	fp = fopen(file.Qobs, "r");
+    if (file.qobs) {
+	if ((fp = fopen(file.qobs, "r")) == NULL)
+	    G_fatal_error(_("Unable to open input file <%s>"), file.qobs);
 
-	misc.Qobs = (double *)G_malloc(input.ntimestep * sizeof(double));
+	misc.Qobs = (double *)G_malloc(input.ntimesteps * sizeof(double));
 
-	for (i = 0; i < input.ntimestep && !feof(fp);) {
+	for (i = 0; i < input.ntimesteps && !feof(fp);) {
 	    get_line(fp, buf);
 
 	    if (sscanf(buf, "%lf", &(misc.Qobs[i])) == 1)
 		i++;
 	}
 
-	input.ntimestep = (input.ntimestep < i ? input.ntimestep : i);
+	input.ntimesteps = (input.ntimesteps < i ? input.ntimesteps : i);
 	fclose(fp);
     }
 
 
-    if (!(misc.timestep > 0 && misc.timestep < input.ntimestep + 1))
+    if (!(misc.timestep > 0 && misc.timestep < input.ntimesteps + 1))
 	misc.timestep = 0;
-    if (!(misc.idxclass > 0 && misc.idxclass < misc.nidxclass + 1))
-	misc.idxclass = 0;
-
+    if (!(misc.topidxclass > 0 && misc.topidxclass < misc.ntopidxclasses + 1))
+	misc.topidxclass = 0;
 
     return;
 }
 
-
-void write_outputs(void)
+void write_output(void)
 {
     FILE *fp;
     time_t tloc;
@@ -203,14 +215,15 @@ void write_outputs(void)
     ltime->tm_mon++;
 
 
-    fp = fopen(file.output, "w");
+    if ((fp = fopen(file.output, "w")) == NULL)
+	G_fatal_error(_("Unable to open output file <%s>"), file.output);
 
     fprintf(fp, "# r.topmodel output file for \"%s\"\n", params.name);
     fprintf(fp, "# Run time: %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
 	    ltime->tm_year, ltime->tm_mon, ltime->tm_mday,
 	    ltime->tm_hour, ltime->tm_min, ltime->tm_sec);
     fprintf(fp, "#\n");
-    if (file.Qobs) {
+    if (file.qobs) {
 	fprintf(fp, "# %-15s Model efficiency\n", "Em:");
 	fprintf(fp, "# %-15s Peak observed Q\n"
 		"# %77s\n", "Qobs_peak:", "[m^3/timestep]");
@@ -225,13 +238,13 @@ void write_outputs(void)
 	    "# %77s\n", "tt_peak:", "[timestep]");
     fprintf(fp, "# %-15s Mean simulated Q\n"
 	    "# %77s\n", "Qt_mean:", "[m^3/timestep]");
-    fprintf(fp, "# %-15s Number of non-NULL cells\n", "ncell:");
+    fprintf(fp, "# %-15s Number of non-NULL cells\n", "ncells:");
     fprintf(fp, "# %-15s Number of topographic index classes\n",
-	    "nidxclass:");
+	    "ntopidxclasses:");
     fprintf(fp, "# %-15s Number of delay timesteps (delay time between "
-	    "rainfall and\n#\t\t\tflow response)\n", "ndelay:");
+	    "rainfall and\n#\t\t\tflow response)\n", "ndelays:");
     fprintf(fp, "# %-15s Number of reach timesteps "
-	    "(time of concentration)\n", "nreach:");
+	    "(time of concentration)\n", "nreaches:");
     fprintf(fp, "# %-15s Areal average of ln(T0) = ln(Te)\n"
 	    "# %77s\n", "lnTe:", "[ln(m^2/timestep)]");
     fprintf(fp, "# %-15s Main channel routing velocity\n"
@@ -267,7 +280,7 @@ void write_outputs(void)
 		"# %77s\n", "fex:", "[m/timestep]");
     }
 
-    if (misc.timestep || misc.idxclass) {
+    if (misc.timestep || misc.topidxclass) {
 	fprintf(fp, "#\n");
 	fprintf(fp, "# %-15s Root zone storage deficit\n"
 		"# %77s\n", "Srz:", "[m]");
@@ -283,59 +296,59 @@ void write_outputs(void)
 
     fprintf(fp, "\n");
 
-    if (file.Qobs) {
-	fprintf(fp, "%-10s ", "Em:");
+    if (file.qobs) {
+	fprintf(fp, "%-16s ", "Em:");
 	if (!Rast_is_d_null_value(&misc.Em))
 	    fprintf(fp, "%10.5lf\n", misc.Em);
 	else
 	    fprintf(fp, "Not resolved due to constant observed Q\n");
-	fprintf(fp, "%-10s %10.3le\n", "Qobs_peak:", misc.Qobs_peak);
-	fprintf(fp, "%-10s %10d\n", "tobs_peak:", misc.tobs_peak);
-	fprintf(fp, "%-10s %10.3le\n", "Qobs_mean:", misc.Qobs_mean);
+	fprintf(fp, "%-16s %10.3le\n", "Qobs_peak:", misc.Qobs_peak);
+	fprintf(fp, "%-16s %10d\n", "tobs_peak:", misc.tobs_peak);
+	fprintf(fp, "%-16s %10.3le\n", "Qobs_mean:", misc.Qobs_mean);
     }
-    fprintf(fp, "%-10s %10.3le\n", "Qt_peak:", misc.Qt_peak);
-    fprintf(fp, "%-10s %10d\n", "tt_peak:", misc.tt_peak);
-    fprintf(fp, "%-10s %10.3le\n", "Qt_mean:", misc.Qt_mean);
-    fprintf(fp, "%-10s %10d\n", "ncell:", misc.ncell);
-    fprintf(fp, "%-10s %10d\n", "nidxclass:", misc.nidxclass);
-    fprintf(fp, "%-10s %10d\n", "ndelay:", misc.ndelay);
-    fprintf(fp, "%-10s %10d\n", "nreach:", misc.nreach);
-    fprintf(fp, "%-10s %10.3le\n", "lnTe:", misc.lnTe);
-    fprintf(fp, "%-10s %10.3le\n", "vch:", misc.vch);
-    fprintf(fp, "%-10s %10.3le\n", "vr:", misc.vr);
-    fprintf(fp, "%-10s %10.3le\n", "lambda:", misc.lambda);
-    fprintf(fp, "%-10s %10.3le\n", "qss:", misc.qss);
-    fprintf(fp, "%-10s %10.3le\n", "qs0:", misc.qs0);
+    fprintf(fp, "%-16s %10.3le\n", "Qt_peak:", misc.Qt_peak);
+    fprintf(fp, "%-16s %10d\n", "tt_peak:", misc.tt_peak);
+    fprintf(fp, "%-16s %10.3le\n", "Qt_mean:", misc.Qt_mean);
+    fprintf(fp, "%-16s %10d\n", "ncells:", misc.ncells);
+    fprintf(fp, "%-16s %10d\n", "ntopidxclasses:", misc.ntopidxclasses);
+    fprintf(fp, "%-16s %10d\n", "ndelays:", misc.ndelays);
+    fprintf(fp, "%-16s %10d\n", "nreaches:", misc.nreaches);
+    fprintf(fp, "%-16s %10.3le\n", "lnTe:", misc.lnTe);
+    fprintf(fp, "%-16s %10.3le\n", "vch:", misc.vch);
+    fprintf(fp, "%-16s %10.3le\n", "vr:", misc.vr);
+    fprintf(fp, "%-16s %10.3le\n", "lambda:", misc.lambda);
+    fprintf(fp, "%-16s %10.3le\n", "qss:", misc.qss);
+    fprintf(fp, "%-16s %10.3le\n", "qs0:", misc.qs0);
     fprintf(fp, "\n");
 
-
     fprintf(fp, "%10s\n", "tch");
     for (i = 0; i < params.nch; i++)
 	fprintf(fp, "%10.3le\n", misc.tch[i]);
+    fprintf(fp, "\n");
 
     fprintf(fp, "%10s\n", "Ad");
-    for (i = 0; i < misc.nreach; i++)
+    for (i = 0; i < misc.nreaches; i++)
 	fprintf(fp, "%10.3le\n", misc.Ad[i]);
-
+    fprintf(fp, "\n");
 
     st = et = si = ei = 0;
-    if (misc.timestep || misc.idxclass) {
+    if (misc.timestep || misc.topidxclass) {
 	if (misc.timestep) {
 	    st = misc.timestep - 1;
 	    et = misc.timestep;
 	}
 	else {
 	    st = 0;
-	    et = input.ntimestep;
+	    et = input.ntimesteps;
 	}
 
-	if (misc.idxclass) {
-	    si = misc.idxclass - 1;
-	    ei = misc.idxclass;
+	if (misc.topidxclass) {
+	    si = misc.topidxclass - 1;
+	    ei = misc.topidxclass;
 	}
 	else {
 	    si = 0;
-	    ei = misc.nidxclass;
+	    ei = misc.ntopidxclasses;
 	}
     }
 
@@ -345,31 +358,32 @@ void write_outputs(void)
 	fprintf(fp, " %10s %10s", "f", "fex");
     fprintf(fp, "\n");
 
-    for (i = 0; i < input.ntimestep; i++) {
+    for (i = 0; i < input.ntimesteps; i++) {
 	fprintf(fp, "%10d %10.3le %10.3le %10.3le %10.3le %10.3le "
 		"%10.3le", i + 1, misc.Qt[i],
-		misc.qt[i][misc.nidxclass],
-		misc.qo[i][misc.nidxclass], misc.qs[i],
-		misc.qv[i][misc.nidxclass], misc.S_mean[i]);
+		misc.qt[i][misc.ntopidxclasses],
+		misc.qo[i][misc.ntopidxclasses], misc.qs[i],
+		misc.qv[i][misc.ntopidxclasses], misc.S_mean[i]);
 	if (params.infex)
 	    fprintf(fp, " %10.3le %10.3le", misc.f[i], misc.fex[i]);
 	fprintf(fp, "\n");
     }
 
-    if (misc.timestep || misc.idxclass) {
+    if (misc.timestep || misc.topidxclass) {
+	fprintf(fp, "\n");
 	fprintf(fp, "Given ");
 	if (misc.timestep)
 	    fprintf(fp, "timestep: %5d", misc.timestep);
-	if (misc.timestep && misc.idxclass)
+	if (misc.timestep && misc.topidxclass)
 	    fprintf(fp, ", ");
-	if (misc.idxclass)
-	    fprintf(fp, "idxclass: %5d", misc.idxclass);
+	if (misc.topidxclass)
+	    fprintf(fp, "topidxclass: %5d", misc.topidxclass);
 	fprintf(fp, "\n");
 
-	if (misc.timestep && !misc.idxclass) {
-	    fprintf(fp, "%10s ", "idxclass");
+	if (misc.timestep && !misc.topidxclass) {
+	    fprintf(fp, "%10s ", "topidxclass");
 	}
-	else if (misc.idxclass && !misc.timestep) {
+	else if (misc.topidxclass && !misc.timestep) {
 	    fprintf(fp, "%10s ", "timestep");
 	}
 
@@ -378,10 +392,10 @@ void write_outputs(void)
 
 	for (i = st; i < et; i++)
 	    for (j = si; j < ei; j++) {
-		if (misc.timestep && !misc.idxclass) {
+		if (misc.timestep && !misc.topidxclass) {
 		    fprintf(fp, "%10d ", j + 1);
 		}
-		else if (misc.idxclass && !misc.timestep) {
+		else if (misc.topidxclass && !misc.timestep) {
 		    fprintf(fp, "%10d ", i + 1);
 		}
 
@@ -397,6 +411,5 @@ void write_outputs(void)
 
     fclose(fp);
 
-
     return;
 }

raster/r.topmodel/global.h

@@ -1,11 +1,4 @@
 #include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <unistd.h>
-#include <math.h>
-#include <time.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
 
 #define	FILL		0x1
 #define	DIR		0x2
@@ -21,90 +14,98 @@
 #define	NTERMS		10
 
 
-/* check_ready.c */
-int check_ready(void);
-int check_required(void);
-int check_names(void);
-int check_io(void);
-
-/* misc.c */
-void gregion(void);
-void depressionless(void);
-void basin_elevation(void);
-void top_index(void);
-
 /* file_io.c */
 void get_line(FILE * fp, char *buffer);
-void read_inputs(void);
-void write_outputs(void);
+void read_input(void);
+void write_output(void);
 
 /* topmodel.c */
-double get_lambda(void);
+void create_topidxstats(char *topidx, int ntopidxclasses, char *outtopidxstats);
+double calculate_lambda(void);
 void initialize(void);
-void implement(void);
-double get_Em(void);
-void others(void);
-void topmodel(void);
+void calculate_flows(void);
+double calculate_Em(void);
+void calculate_others(void);
+void run_topmodel(void);
 
 /* infiltration.c */
-double get_f(double t, double R);
+double calculate_f(double t, double R);
 
 
 /* Topographic index statistics file */
-struct idxstats
+struct topidxstats
 {
-    /* misc.nidxclass's */
-    double *atb, *Aatb_r;
+    /* misc.ntopidxclasses */
+    double *atb;
+    double *Aatb_r;
 };
 
 /* Parameters file */
 struct params
 {
     char *name;
-    double A, qs0, lnTe, m, Sr0, Srmax, td, vch, vr;
+    double A;
+    double qs0;
+    double lnTe;
+    double m;
+    double Sr0;
+    double Srmax;
+    double td;
+    double vch;
+    double vr;
     int infex;
-    double K0, psi, dtheta;
+    double K0;
+    double psi;
+    double dtheta;
     int nch;
     /* params.nch's */
-    double *d, *Ad_r;
+    double *d;
+    double *Ad_r;
 };
 
 /* Input file */
 struct input
 {
-    int ntimestep;
+    int ntimesteps;
     double dt;
     /* input.ntimestep's */
-    double *R, *Ep;
-};
-
-/* Map names */
-struct map
-{
-    char *basin, *elev, *belev, *fill, *dir, *topidx;
+    double *R;
+    double *Ep;
 };
 
 /* File names */
 struct file
 {
-    char *idxstats, *params, *input, *output, *Qobs;
+    char *params;
+    char *topidxstats;
+    char *input;
+    char *output;
+    char *qobs;
 };
 
 /* Miscellaneous TOPMODEL variables */
 struct misc
 {
     /* Number of non-null cells */
-    int ncell;
+    int ncells;
     /* Number of topographic index classes */
-    int nidxclass;
+    int ntopidxclasses;
     /* Model efficiency */
     double Em;
-    int ndelay, nreach;
-    double lnTe, vch, vr;
+    int ndelays;
+    int nreaches;
+    double lnTe;
+    double vch;
+    double vr;
     double lambda;
-    double qss, qs0;
-    double Qobs_peak, Qt_peak, Qobs_mean, Qt_mean;
-    int tobs_peak, tt_peak;
+    double qss;
+    double qs0;
+    double Qobs_peak;
+    double Qt_peak;
+    double Qobs_mean;
+    double Qt_mean;
+    int tobs_peak;
+    int tt_peak;
     /* params.nch's */
     double *tch;
     /* misc.nreach's */
@@ -116,35 +117,32 @@ struct misc
     double *S_mean;
     double *f;
     double *fex;
-    /* input.ntimestep * (misc.nidxclass + 1)'s */
-    double **qt, **qo, **qv;
-    /* input.ntimestep * misc.nidxclass's */
-    double **Srz, **Suz;
+    /* input.ntimestep * (misc.ntopidxclasses + 1)'s */
+    double **qt;
+    double **qo;
+    double **qv;
+    /* input.ntimestep * misc.ntopidxclassess */
+    double **Srz;
+    double **Suz;
     double **S;
     double **Ea;
     double **ex;
     /* Miscellaneous variables */
-    int timestep, idxclass;
+    int timestep;
+    int topidxclass;
 };
 
-struct flg
-{
-    /* Input flag */
-    char input;
-    /* Overwrite flag */
-    char overwr;
-    /* Overwrite list */
-    int overwrlist;
-};
+#ifdef _MAIN_C_
+#define GLOBAL
+#else
+#define GLOBAL extern
+#endif
 
-extern struct idxstats idxstats;
-extern struct params params;
-extern struct input input;
-extern struct map map;
-extern struct file file;
-extern struct misc misc;
-extern struct flg flg;
+GLOBAL struct params params;
+GLOBAL struct topidxstats topidxstats;
+GLOBAL struct input input;
+GLOBAL struct file file;
+GLOBAL struct misc misc;
 
 /* Miscellaneous variables */
-extern const char *mapset;
-extern char buf[BUFSIZE];
+GLOBAL char buf[BUFSIZE];

raster/r.topmodel/infiltration.c

@@ -1,8 +1,9 @@
+#include <math.h>
+#include <grass/raster.h>
 #include "global.h"
 
-
 /* The Green-and-Ampt Model */
-double get_f(double t, double R)
+double calculate_f(double t, double R)
 {
     static double cumf = 0.0, f_ = 0.0;
     static char ponding = 0;

raster/r.topmodel/main.c

@@ -4,7 +4,7 @@
  * TMOD9502.FOR Author: Keith Beven <k.beven@lancaster.ac.uk>
  *                      http://www.es.lancs.ac.uk/hfdg/topmodel.html
  *
- *      Copyright (C) 2000, 2010 by the GRASS Development Team
+ *      Copyright (C) 2000, 2010, 2013 by the GRASS Development Team
  *      Author: Huidae Cho <grass4u@gmail.com>
  *              Hydro Laboratory, Kyungpook National University
  *              South Korea
@@ -14,46 +14,29 @@
  *
  */
 
+#define _MAIN_C_
 #include <stdio.h>
 #include <stdlib.h>
+#include <grass/gis.h>
 #include <grass/glocale.h>
 #include "global.h"
 
-struct idxstats idxstats;
-struct params params;
-struct input input;
-struct map map;
-struct file file;
-struct misc misc;
-struct flg flg;
-const char *mapset;
-
 int main(int argc, char **argv)
 {
-
     struct GModule *module;
     struct
     {
-	struct Option *basin;
-	struct Option *elev;
-	struct Option *fill;
-	struct Option *dir;
-	struct Option *belev;
-	struct Option *topidx;
-	struct Option *nidxclass;
-	struct Option *idxstats;
 	struct Option *params;
+	struct Option *topidxstats;
 	struct Option *input;
 	struct Option *output;
-	struct Option *Qobs;
+	struct Option *qobs;
 	struct Option *timestep;
-	struct Option *idxclass;
-    } param;
-
-    struct
-    {
-	struct Flag *input;
-    } flag;
+	struct Option *topidxclass;
+	struct Option *topidx;
+	struct Option *ntopidxclasses;
+	struct Option *outtopidxstats;
+    } params;
 
     /* Initialize GRASS and parse command line */
     G_gisinit(argv[0]);
@@ -65,154 +48,115 @@ int main(int argc, char **argv)
 	_("Simulates TOPMODEL which is a physically based hydrologic model.");
 
     /* Parameter definitions */
-    param.elev = G_define_standard_option(G_OPT_R_ELEV);
-    param.elev->required = NO;
-    param.elev->guisection = _("Input");
-
-    param.basin = G_define_standard_option(G_OPT_R_INPUT);
-    param.basin->key = "basin";
-    param.basin->label =
-	_("Name of input basin raster map");
-    param.basin->description = _("Created by r.water.outlet (MASK)");
-    param.basin->required = NO;
-    param.basin->guisection = _("Input");
-
-    param.fill = G_define_standard_option(G_OPT_R_OUTPUT);
-    param.fill->key = "depressionless";
-    param.fill->description = _("Name for output depressionless elevation raster map");
-    param.fill->required = NO;
-    param.fill->guisection = _("Output");
-
-    param.dir = G_define_standard_option(G_OPT_R_OUTPUT);
-    param.dir->key = "direction";
-    param.dir->description =
-	_("Name for output flow direction map for depressionless elevation raster map");
-    param.dir->required = NO;
-    param.dir->guisection = _("Output");
-
-    param.belev = G_define_standard_option(G_OPT_R_OUTPUT);
-    param.belev->key = "basin_elevation";
-    param.belev->label = _("Name for output basin elevation raster map (o/i)");
-    param.belev->description = _("MASK applied");
-    param.belev->required = NO;
-    param.belev->guisection = _("Output");
-
-    param.topidx = G_define_standard_option(G_OPT_R_OUTPUT);
-    param.topidx->key = "topidx";
-    param.topidx->label =
-	_("Name for output topographic index ln(a/tanB) raster map");
-    param.topidx->description = _("MASK applied");
-    param.topidx->required = NO;
-    param.topidx->guisection = _("Output");
-
-    param.nidxclass = G_define_option();
-    param.nidxclass->key = "nidxclass";
-    param.nidxclass->description =
-	_("Number of topographic index classes");
-    param.nidxclass->type = TYPE_INTEGER;
-    param.nidxclass->required = NO;
-    param.nidxclass->answer = "30";
-    param.nidxclass->guisection = _("Parameters");
-    
-    param.idxstats = G_define_standard_option(G_OPT_F_INPUT);
-    param.idxstats->key = "idxstats";
-    param.idxstats->description =
-	_("Name of topographic index statistics file (o/i)");
-    
-    param.params = G_define_standard_option(G_OPT_F_INPUT);
-    param.params->key = "parameters";
-    param.params->description = _("Name of TOPMODEL parameters file");
-    
-    param.input = G_define_standard_option(G_OPT_F_INPUT);
-    param.input->description =
+    params.params = G_define_standard_option(G_OPT_F_INPUT);
+    params.params->key = "parameters";
+    params.params->description = _("Name of TOPMODEL parameters file");
+
+    params.topidxstats = G_define_standard_option(G_OPT_F_INPUT);
+    params.topidxstats->key = "topidxstats";
+    params.topidxstats->description =
+	_("Name of topographic index statistics file");
+
+    params.input = G_define_standard_option(G_OPT_F_INPUT);
+    params.input->description =
 	_("Name of rainfall and potential evapotranspiration data file");
 
-    param.output = G_define_standard_option(G_OPT_F_OUTPUT);
-    param.output->description = _("Name for output file");
-
-    param.Qobs = G_define_standard_option(G_OPT_F_OUTPUT);
-    param.Qobs->key = "qobs";
-    param.Qobs->description = _("Name for observed flow file");
-    param.Qobs->required = NO;
-    param.Qobs->guisection = _("Output");
-    
-    param.timestep = G_define_option();
-    param.timestep->key = "timestep";
-    param.timestep->description =
-	_("Time step");
-    param.timestep->type = TYPE_INTEGER;
-    param.timestep->required = NO;
-    param.timestep->guisection = _("Parameters");
-    
-    param.idxclass = G_define_option();
-    param.idxclass->key = "idxclass";
-    param.idxclass->description =
-	_("Topographic index class");
-    param.idxclass->type = TYPE_INTEGER;
-    param.idxclass->required = NO;
-    param.idxclass->guisection = _("Parameters");
-
-    /* Flag definitions */
-    flag.input = G_define_flag();
-    flag.input->key = 'i';
-    flag.input->description = _("Input data given for (o/i)");
+    params.output = G_define_standard_option(G_OPT_F_OUTPUT);
+    params.output->description = _("Name for output file");
+
+    params.qobs = G_define_standard_option(G_OPT_F_INPUT);
+    params.qobs->key = "qobs";
+    params.qobs->description = _("Name of observed flow file");
+    params.qobs->required = NO;
+
+    params.timestep = G_define_option();
+    params.timestep->key = "timestep";
+    params.timestep->label = _("Time step");
+    params.timestep->description = _("Generate output for this time step.");
+    params.timestep->type = TYPE_INTEGER;
+    params.timestep->required = NO;
+
+    params.topidxclass = G_define_option();
+    params.topidxclass->key = "topidxclass";
+    params.topidxclass->label = _("Topographic index class");
+    params.topidxclass->description =
+	_("Generate output for this topographic index class.");
+    params.topidxclass->type = TYPE_INTEGER;
+    params.topidxclass->required = NO;
+
+    params.topidx = G_define_standard_option(G_OPT_R_INPUT);
+    params.topidx->key = "topidx";
+    params.topidx->label =
+	_("Name of input topographic index ln(a/tanB) raster map");
+    params.topidx->description =
+	    _("Must be clipped to the catchment boundary. Used for generating outtopidxstats.");
+    params.topidx->required = NO;
+    params.topidx->guisection = _("Preprocess");
+
+    params.ntopidxclasses = G_define_option();
+    params.ntopidxclasses->key = "ntopidxclasses";
+    params.ntopidxclasses->label = _("Number of topographic index classes");
+    params.ntopidxclasses->description =
+	    _("Used for generating outtopidxstats.");
+    params.ntopidxclasses->type = TYPE_INTEGER;
+    params.ntopidxclasses->required = NO;
+    params.ntopidxclasses->answer = "30";
+    params.ntopidxclasses->guisection = _("Preprocess");
+
+    params.outtopidxstats = G_define_standard_option(G_OPT_F_OUTPUT);
+    params.outtopidxstats->key = "outtopidxstats";
+    params.outtopidxstats->label =
+	_("Name for output topographic index statistics file");
+    params.outtopidxstats->description =
+	    _("Requires topidx and ntopidxclasses.");
+    params.outtopidxstats->required = NO;
+    params.outtopidxstats->guisection = _("Preprocess");
 
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-    /* Store given parameters and flags */
-    map.basin = param.basin->answer;
-    map.elev = param.elev->answer;
-    map.belev = param.belev->answer;
-    map.fill = param.fill->answer;
-    map.dir = param.dir->answer;
-    map.topidx = param.topidx->answer;
-    file.idxstats = param.idxstats->answer;
-    file.params = param.params->answer;
-    file.input = param.input->answer;
-    file.output = param.output->answer;
-    file.Qobs = param.Qobs->answer;
-
-    misc.nidxclass = atoi(param.nidxclass->answer);
-
-    if (!param.timestep->answer)
-	param.timestep->answer = "0";
-    if (!param.idxclass->answer)
-	param.idxclass->answer = "0";
-
-    misc.timestep = atoi(param.timestep->answer);
-    misc.idxclass = atoi(param.idxclass->answer);
-
-    flg.input = flag.input->answer;
-    flg.overwr = module->overwrite;
-
-    mapset = G_mapset();
-
-    /* Check run conditions */
-    if (check_ready())
-	exit(1);
-
-    /* Adjust cell header */
-    gregion();
-
-    /* Create required maps */
-    if (!flg.input) {
-	if (map.fill)
-	    depressionless();
-	basin_elevation();
+    /* Store given parameters */
+    file.params = params.params->answer;
+    file.topidxstats = params.topidxstats->answer;
+    file.input = params.input->answer;
+    file.output = params.output->answer;
+    file.qobs = params.qobs->answer;
+
+    if (!params.timestep->answer)
+	params.timestep->answer = "0";
+    misc.timestep = atoi(params.timestep->answer);
+
+    if (!params.topidxclass->answer)
+	params.topidxclass->answer = "0";
+    misc.topidxclass = atoi(params.topidxclass->answer);
+
+    if (params.topidx->answer && params.outtopidxstats->answer) {
+	char *topidx;
+	int ntopidxclasses;
+	char *outtopidxstats;
+	
+	topidx = params.topidx->answer;
+	ntopidxclasses = atoi(params.ntopidxclasses->answer);
+	outtopidxstats = params.outtopidxstats->answer;
+
+	if (ntopidxclasses <= 0)
+	    G_fatal_error(_("ntopidxclasses must be positive."));
+
+	create_topidxstats(topidx, ntopidxclasses, outtopidxstats);
+    } else if (params.topidx->answer) {
+	G_warning(_("Ignoring topidx because outtopidxstats is not specified."));
+    } else if (params.outtopidxstats->answer) {
+	G_warning(_("Ignoring outtopidxstats because topidx is not specified."));
     }
 
-    top_index();
-
     /* Read required files */
-    read_inputs();
+    read_input();
 
     /* Implement TOPMODEL */
-    topmodel();
+    run_topmodel();
 
     /* Write outputs */
-    write_outputs();
-
+    write_output();
 
     exit(0);
 }

raster/r.topmodel/r.topmodel.html

@@ -2,37 +2,43 @@
 
 <b><em>r.topmodel</em></b> simulates TOPMODEL which is a physically based
 hydrologic model.
-<p>Note: (i) means input; (o) means output; (o/i) means input or output
-<p>The <b>-i</b> flag indicates that input data are given for (o/i).  Without this
-flag, all inputs (i) and intermediate outputs (o/i) should be given.  For
-example, [belevation] map will be created from [elevation] and [basin] in every
-run.  However, given the same [elevation] and [basin], [belevation] output will
-be the same all the time, so r.topmodel can directly take [belevation] as an
-input with this flag to save time.
-<p>
+
 <h3>Selected Parameters:</h3>
 
 <dl>
-<dt><b>depressionless</b> map is created as follows:</dt>
-<dd><div class="code"><pre>
-r.fill.dir input=elevation elev=depressionless dir=direction type=grass
-</pre></div>
-This option can be omitted if [elevation] map is already depressionless.
+<dt><b>qobs</b></dt>
+<dd>Compare simulated flows with observed flows and calculate the model
+efficiency.
 </dd>
 
-<dt><b>belevation</b> map is created from [elevation] with [basin] mask applied:</dt>
-<dd><div class="code"><pre>
-r.mapcalc "belevation = if(basin == 0 || isnull(basin), null(), elevation)"
-</pre></div></dd>
+<dt><b>timestep</b></dt>
+<dd>
+If a time step is specified, output will be generated for the specific time
+step in addition to the summary and total flows at the outlet. This parameter
+can be combined with topidxclass to specify a time step and topographic index
+class at the same time. If no topidxclass is given, output will be generated
+for all the topographic index classes.
+</dd>
 
-<dt><b>topidx</b> map is created as follows:</dt>
-<dd><div class="code"><pre>
-r.topidx input=elevation output=topidx
-</pre></div></dd>
+<dt><b>topidxclass</b></dt>
+<dd>
+If a topographic index class is specified, output will be generated for the
+given topographic index class. This parameter can be combined with timestep. If
+no timestep is given, output will be generated for all the time steps.
+</dd>
 
-<dt><b>qobs</b></dt>
-<dd>Compare simulated flows with observed flows and calculate model
-efficiency.
+<dt><b>topidx</b>, <b>ntopidxclasses</b>, <b>outtopidxstats</b></dt>
+<dd>
+The <b>topidx</b> map can optionally be used for creating a new topographic
+index statistics file. This map has to be already clipped to the catchment
+boundary. The entire range of topographic index values will be divided into
+<b>ntopidxclasses</b> and the number of cells in each class will be reported in
+the <b>outtopidxstats</b> file using the following command:
+<div class="code"><pre>
+r.stats -Anc input=[topidx] output=[outtopidxstats] nsteps=[ntopidxclasses]
+</pre></div>
+These three parameters can be omitted unless new topidxstats file needs to be
+created.
 </dd>
 </dl>
 
@@ -42,7 +48,8 @@ efficiency.
 K. Beven, R. Lamb, P. Quinn, R. Romanowicz, and J. Freer.
 TOPMODEL, in V.P. Singh (Ed.). Computer Models of Watershed Hydrology.
 Water Resources Publications, 1995.
-<p>S.C. Liaw, Streamflow simulation using a physically based hydrologic
+<p>
+S.C. Liaw, Streamflow simulation using a physically based hydrologic
 model in humid forested watersheds (Dissertation).
 Colorado State University, CO. p163, 1988.
 
@@ -58,10 +65,9 @@ Colorado State University, CO. p163, 1988.
 
 <h2>AUTHORS</h2>
 
-Main algorithm sources are rewritten in C based on TMOD9502.FOR.
-<br>
-Thanks to Keith Beven.
-<p>GRASS port by <a href="mailto:grass4u@gmail com">Huidae Cho</a><br>
+<a href="mailto:grass4u@gmail com">Huidae Cho</a><br>
 Hydro Laboratory, Kyungpook National University, South Korea
+<p>
+Based on TMOD9502.FOR by Keith Beven.
 
 <p><i>Last changed: $Date$</i>

raster/r.topmodel/topmodel.c

@@ -1,31 +1,48 @@
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
+#include <grass/spawn.h>
 #include "global.h"
 
+void create_topidxstats(char *topidx, int ntopidxclasses, char *outtopidxstats)
+{
+    char input[GPATH_MAX];
+    char output[GPATH_MAX];
+    char nsteps[32];
+
+    sprintf(input, "input=%s", topidx);
+    sprintf(output, "output=%s", outtopidxstats);
+    sprintf(nsteps, "nsteps=%d", ntopidxclasses);
 
-double get_lambda(void)
+    G_message("Creating topographic index statistics file...");
+    G_verbose_message("r.stats -Anc %s %s %s ...", input, output, nsteps);
+
+    if (G_spawn("r.stats", "r.stats", "-Anc", input, output, nsteps, NULL) != 0)
+	G_fatal_error("r.stats failed");
+}
+
+double calculate_lambda(void)
 {
     int i;
     double retval;
 
 
     retval = 0.0;
-    for (i = 1; i < misc.nidxclass; i++)
-	retval += idxstats.Aatb_r[i] *
-	    (idxstats.atb[i] + idxstats.atb[i - 1]) / 2.0;
+    for (i = 1; i < misc.ntopidxclasses; i++)
+	retval += topidxstats.Aatb_r[i] *
+	    (topidxstats.atb[i] + topidxstats.atb[i - 1]) / 2.0;
 
 
     return retval;
 }
 
-
 void initialize(void)
 {
     int i, j, t;
     double A1, A2;
 
 
-    misc.lambda = get_lambda();
+    misc.lambda = calculate_lambda();
     misc.lnTe = params.lnTe + log(input.dt);
     misc.vch = params.vch * input.dt;
     misc.vr = params.vr * input.dt;
@@ -37,16 +54,16 @@ void initialize(void)
     for (i = 1; i < params.nch; i++)
 	misc.tch[i] = misc.tch[0] + (params.d[i] - params.d[0]) / misc.vr;
 
-    misc.nreach = (int)misc.tch[params.nch - 1];
-    if ((double)misc.nreach < misc.tch[params.nch - 1])
-	misc.nreach++;
-    misc.ndelay = (int)misc.tch[0];
+    misc.nreaches = (int)misc.tch[params.nch - 1];
+    if ((double)misc.nreaches < misc.tch[params.nch - 1])
+	misc.nreaches++;
+    misc.ndelays = (int)misc.tch[0];
 
-    misc.nreach -= misc.ndelay;
+    misc.nreaches -= misc.ndelays;
 
-    misc.Ad = (double *)G_malloc(misc.nreach * sizeof(double));
-    for (i = 0; i < misc.nreach; i++) {
-	t = misc.ndelay + i + 1;
+    misc.Ad = (double *)G_malloc(misc.nreaches * sizeof(double));
+    for (i = 0; i < misc.nreaches; i++) {
+	t = misc.ndelays + i + 1;
 	if (t > misc.tch[params.nch - 1]) {
 	    misc.Ad[i] = 1.0;
 	}
@@ -66,47 +83,46 @@ void initialize(void)
 
     A1 = misc.Ad[0];
     misc.Ad[0] *= params.A;
-    for (i = 1; i < misc.nreach; i++) {
+    for (i = 1; i < misc.nreaches; i++) {
 	A2 = misc.Ad[i];
 	misc.Ad[i] = A2 - A1;
 	A1 = A2;
 	misc.Ad[i] *= params.A;
     }
 
-    misc.Srz = (double **)G_malloc(input.ntimestep * sizeof(double *));
-    misc.Suz = (double **)G_malloc(input.ntimestep * sizeof(double *));
-    for (i = 0; i < input.ntimestep; i++) {
-	misc.Srz[i] = (double *)G_malloc(misc.nidxclass * sizeof(double));
-	misc.Suz[i] = (double *)G_malloc(misc.nidxclass * sizeof(double));
+    misc.Srz = (double **)G_malloc(input.ntimesteps * sizeof(double *));
+    misc.Suz = (double **)G_malloc(input.ntimesteps * sizeof(double *));
+    for (i = 0; i < input.ntimesteps; i++) {
+	misc.Srz[i] = (double *)G_malloc(misc.ntopidxclasses * sizeof(double));
+	misc.Suz[i] = (double *)G_malloc(misc.ntopidxclasses * sizeof(double));
     }
 
-    for (i = 0; i < misc.nidxclass; i++) {
+    for (i = 0; i < misc.ntopidxclasses; i++) {
 	misc.Srz[0][i] = params.Sr0;
 	misc.Suz[0][i] = 0.0;
     }
 
-    misc.S_mean = (double *)G_malloc(input.ntimestep * sizeof(double));
+    misc.S_mean = (double *)G_malloc(input.ntimesteps * sizeof(double));
     misc.S_mean[0] = -params.m * log(misc.qs0 / misc.qss);
 
-    misc.Qt = (double *)G_malloc(input.ntimestep * sizeof(double));
-    for (i = 0; i < input.ntimestep; i++)
+    misc.Qt = (double *)G_malloc(input.ntimesteps * sizeof(double));
+    for (i = 0; i < input.ntimesteps; i++)
 	misc.Qt[i] = 0.0;
 
-    for (i = 0; i < misc.ndelay; i++)
+    for (i = 0; i < misc.ndelays; i++)
 	misc.Qt[i] = misc.qs0 * params.A;
 
     A1 = 0.0;
-    for (i = 0; i < misc.nreach; i++) {
+    for (i = 0; i < misc.nreaches; i++) {
 	A1 += misc.Ad[i];
-	misc.Qt[misc.ndelay + i] = misc.qs0 * (params.A - A1);
+	misc.Qt[misc.ndelays + i] = misc.qs0 * (params.A - A1);
     }
 
 
     return;
 }
 
-
-void implement(void)
+void calculate_flows(void)
 {
     int i, j, k;
     double Aatb_r;
@@ -114,38 +130,38 @@ void implement(void)
     double _qo, _qv;
 
 
-    misc.S = (double **)G_malloc(input.ntimestep * sizeof(double *));
-    misc.Ea = (double **)G_malloc(input.ntimestep * sizeof(double *));
-    misc.ex = (double **)G_malloc(input.ntimestep * sizeof(double *));
+    misc.S = (double **)G_malloc(input.ntimesteps * sizeof(double *));
+    misc.Ea = (double **)G_malloc(input.ntimesteps * sizeof(double *));
+    misc.ex = (double **)G_malloc(input.ntimesteps * sizeof(double *));
 
-    misc.qt = (double **)G_malloc(input.ntimestep * sizeof(double *));
-    misc.qo = (double **)G_malloc(input.ntimestep * sizeof(double *));
-    misc.qv = (double **)G_malloc(input.ntimestep * sizeof(double *));
+    misc.qt = (double **)G_malloc(input.ntimesteps * sizeof(double *));
+    misc.qo = (double **)G_malloc(input.ntimesteps * sizeof(double *));
+    misc.qv = (double **)G_malloc(input.ntimesteps * sizeof(double *));
 
-    misc.qs = (double *)G_malloc(input.ntimestep * sizeof(double));
-    misc.f = (double *)G_malloc(input.ntimestep * sizeof(double));
-    misc.fex = (double *)G_malloc(input.ntimestep * sizeof(double));
+    misc.qs = (double *)G_malloc(input.ntimesteps * sizeof(double));
+    misc.f = (double *)G_malloc(input.ntimesteps * sizeof(double));
+    misc.fex = (double *)G_malloc(input.ntimesteps * sizeof(double));
 
-    for (i = 0; i < input.ntimestep; i++) {
-	misc.S[i] = (double *)G_malloc(misc.nidxclass * sizeof(double));
-	misc.Ea[i] = (double *)G_malloc(misc.nidxclass * sizeof(double));
-	misc.ex[i] = (double *)G_malloc(misc.nidxclass * sizeof(double));
+    for (i = 0; i < input.ntimesteps; i++) {
+	misc.S[i] = (double *)G_malloc(misc.ntopidxclasses * sizeof(double));
+	misc.Ea[i] = (double *)G_malloc(misc.ntopidxclasses * sizeof(double));
+	misc.ex[i] = (double *)G_malloc(misc.ntopidxclasses * sizeof(double));
 
-	misc.qt[i] = (double *)G_malloc((misc.nidxclass + 1) *
+	misc.qt[i] = (double *)G_malloc((misc.ntopidxclasses + 1) *
 					sizeof(double));
-	misc.qo[i] = (double *)G_malloc((misc.nidxclass + 1) *
+	misc.qo[i] = (double *)G_malloc((misc.ntopidxclasses + 1) *
 					sizeof(double));
-	misc.qv[i] = (double *)G_malloc((misc.nidxclass + 1) *
+	misc.qv[i] = (double *)G_malloc((misc.ntopidxclasses + 1) *
 					sizeof(double));
 
-	misc.qt[i][misc.nidxclass] = 0.0;
-	misc.qo[i][misc.nidxclass] = 0.0;
-	misc.qv[i][misc.nidxclass] = 0.0;
+	misc.qt[i][misc.ntopidxclasses] = 0.0;
+	misc.qo[i][misc.ntopidxclasses] = 0.0;
+	misc.qv[i][misc.ntopidxclasses] = 0.0;
 	misc.qs[i] = 0.0;
 
 	if (params.infex) {
 	    misc.f[i] = input.dt *
-		get_f((i + 1) * input.dt, input.R[i] / input.dt);
+		calculate_f((i + 1) * input.dt, input.R[i] / input.dt);
 	    misc.fex[i] = input.R[i] - misc.f[i];
 	    R = misc.f[i];
 	}
@@ -156,7 +172,7 @@ void implement(void)
 	}
 
 	if (i) {
-	    for (j = 0; j < misc.nidxclass; j++) {
+	    for (j = 0; j < misc.ntopidxclasses; j++) {
 		misc.Srz[i][j] = misc.Srz[i - 1][j];
 		misc.Suz[i][j] = misc.Suz[i - 1][j];
 	    }
@@ -164,13 +180,13 @@ void implement(void)
 
 	misc.qs[i] = misc.qss * exp(-misc.S_mean[i] / params.m);
 
-	for (j = 0; j < misc.nidxclass; j++) {
-	    Aatb_r = (idxstats.Aatb_r[j] +
-		      (j < misc.nidxclass - 1 ? idxstats.Aatb_r[j + 1]
+	for (j = 0; j < misc.ntopidxclasses; j++) {
+	    Aatb_r = (topidxstats.Aatb_r[j] +
+		      (j < misc.ntopidxclasses - 1 ? topidxstats.Aatb_r[j + 1]
 		       : 0.0)) / 2.0;
 
 	    misc.S[i][j] = misc.S_mean[i] +
-		params.m * (misc.lambda - idxstats.atb[j]);
+		params.m * (misc.lambda - topidxstats.atb[j]);
 	    if (misc.S[i][j] < 0.0)
 		misc.S[i][j] = 0.0;
 
@@ -202,7 +218,7 @@ void implement(void)
 		_qv *= Aatb_r;
 	    }
 	    misc.qv[i][j] = _qv;
-	    misc.qv[i][misc.nidxclass] += misc.qv[i][j];
+	    misc.qv[i][misc.ntopidxclasses] += misc.qv[i][j];
 
 	    misc.Ea[i][j] = 0.0;
 	    if (input.Ep[i] > 0.0) {
@@ -216,7 +232,7 @@ void implement(void)
 	    _qo = 0.0;
 	    if (j > 0) {
 		if (misc.ex[i][j] > 0.0)
-		    _qo = idxstats.Aatb_r[j] *
+		    _qo = topidxstats.Aatb_r[j] *
 			(misc.ex[i][j - 1] + misc.ex[i][j]) / 2.0;
 		else if (misc.ex[i][j - 1] > 0.0)
 		    _qo = Aatb_r * misc.ex[i][j - 1] /
@@ -224,24 +240,24 @@ void implement(void)
 			 misc.ex[i][j]) * misc.ex[i][j - 1] / 2.0;
 	    }
 	    misc.qo[i][j] = _qo;
-	    misc.qo[i][misc.nidxclass] += misc.qo[i][j];
+	    misc.qo[i][misc.ntopidxclasses] += misc.qo[i][j];
 
 	    misc.qt[i][j] = misc.qo[i][j] + misc.qs[i];
 	}
-	misc.qo[i][misc.nidxclass] += misc.fex[i];
-	misc.qt[i][misc.nidxclass] = misc.qo[i][misc.nidxclass] + misc.qs[i];
+	misc.qo[i][misc.ntopidxclasses] += misc.fex[i];
+	misc.qt[i][misc.ntopidxclasses] = misc.qo[i][misc.ntopidxclasses] + misc.qs[i];
 
 	misc.S_mean[i] = misc.S_mean[i] +
-	    misc.qs[i] - misc.qv[i][misc.nidxclass];
+	    misc.qs[i] - misc.qv[i][misc.ntopidxclasses];
 
-	if (i + 1 < input.ntimestep)
+	if (i + 1 < input.ntimesteps)
 	    misc.S_mean[i + 1] = misc.S_mean[i];
 
-	for (j = 0; j < misc.nreach; j++) {
-	    k = i + j + misc.ndelay;
-	    if (k > input.ntimestep - 1)
+	for (j = 0; j < misc.nreaches; j++) {
+	    k = i + j + misc.ndelays;
+	    if (k > input.ntimesteps - 1)
 		break;
-	    misc.Qt[k] += misc.qt[i][misc.nidxclass] * misc.Ad[j];
+	    misc.Qt[k] += misc.qt[i][misc.ntopidxclasses] * misc.Ad[j];
 	}
     }
 
@@ -249,9 +265,8 @@ void implement(void)
     return;
 }
 
-
-/* Object function for hydrograph suggested by Servet and Dezetter(1991) */
-double get_Em(void)
+/* Objective function for hydrograph suggested by Servet and Dezetter(1991) */
+double calculate_Em(void)
 {
     int i;
     double Em, numerator, denominator;
@@ -259,14 +274,14 @@ double get_Em(void)
 
     misc.Qobs_mean = 0.0;
     numerator = 0.0;
-    for (i = 0; i < input.ntimestep; i++) {
+    for (i = 0; i < input.ntimesteps; i++) {
 	misc.Qobs_mean += misc.Qobs[i];
 	numerator += pow(misc.Qobs[i] - misc.Qt[i], 2.0);
     }
-    misc.Qobs_mean /= input.ntimestep;
+    misc.Qobs_mean /= input.ntimesteps;
 
     denominator = 0.0;
-    for (i = 0; i < input.ntimestep; i++)
+    for (i = 0; i < input.ntimesteps; i++)
 	denominator += pow(misc.Qobs[i] - misc.Qobs_mean, 2.0);
 
     if (denominator == 0.0) {
@@ -277,29 +292,27 @@ double get_Em(void)
 	Em = 1.0 - numerator / denominator;
     }
 
-
     return Em;
 }
 
-
-void others(void)
+void calculate_others(void)
 {
     int i;
 
 
     misc.Qt_mean = 0.0;
-    for (i = 0; i < input.ntimestep; i++) {
+    for (i = 0; i < input.ntimesteps; i++) {
 	misc.Qt_mean += misc.Qt[i];
 	if (!i || misc.Qt_peak < misc.Qt[i]) {
 	    misc.Qt_peak = misc.Qt[i];
 	    misc.tt_peak = i + 1;
 	}
     }
-    misc.Qt_mean /= input.ntimestep;
+    misc.Qt_mean /= input.ntimesteps;
 
-    if (file.Qobs) {
-	misc.Em = get_Em();
-	for (i = 0; i < input.ntimestep; i++) {
+    if (file.qobs) {
+	misc.Em = calculate_Em();
+	for (i = 0; i < input.ntimesteps; i++) {
 	    if (!i || misc.Qobs_peak < misc.Qobs[i]) {
 		misc.Qobs_peak = misc.Qobs[i];
 		misc.tobs_peak = i + 1;
@@ -307,17 +320,14 @@ void others(void)
 	}
     }
 
-
     return;
 }
 
-
-void topmodel(void)
+void run_topmodel(void)
 {
     initialize();
-    implement();
-    others();
-
+    calculate_flows();
+    calculate_others();
 
     return;
 }