[GRASS-SVN] r58073 - grass/trunk/raster/r.topmodel
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Oct 20 05:19:11 PDT 2013
Author: hcho
Date: 2013-10-20 05:19:10 -0700 (Sun, 20 Oct 2013)
New Revision: 58073
Modified:
grass/trunk/raster/r.topmodel/file_io.c
grass/trunk/raster/r.topmodel/global.h
grass/trunk/raster/r.topmodel/infiltration.c
grass/trunk/raster/r.topmodel/main.c
grass/trunk/raster/r.topmodel/r.topmodel.html
grass/trunk/raster/r.topmodel/topmodel.c
Log:
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.
Modified: grass/trunk/raster/r.topmodel/file_io.c
===================================================================
--- grass/trunk/raster/r.topmodel/file_io.c 2013-10-20 03:57:31 UTC (rev 58072)
+++ grass/trunk/raster/r.topmodel/file_io.c 2013-10-20 12:19:10 UTC (rev 58073)
@@ -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 @@
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 @@
/* 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);
- misc.ncell = 0;
+ topidxstats.atb = NULL;
+ topidxstats.Aatb_r = NULL;
+ misc.ncells = 0;
- for (i = 0; i < misc.nidxclass && !feof(fp);) {
+ for (i = 0, j = 0; !feof(fp);) {
+ double atb;
+ double Aatb_r;
+
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.ntopidxclasses; i++)
+ topidxstats.Aatb_r[i] /= (double)misc.ncells;
- 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++) {
+ 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 @@
/* 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 @@
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 @@
"# %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 @@
"# %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 @@
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 @@
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 @@
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 @@
fclose(fp);
-
return;
}
Modified: grass/trunk/raster/r.topmodel/global.h
===================================================================
--- grass/trunk/raster/r.topmodel/global.h 2013-10-20 03:57:31 UTC (rev 58072)
+++ grass/trunk/raster/r.topmodel/global.h 2013-10-20 12:19:10 UTC (rev 58073)
@@ -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;
+ double *R;
+ double *Ep;
};
-/* Map names */
-struct map
-{
- char *basin, *elev, *belev, *fill, *dir, *topidx;
-};
-
/* 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 @@
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];
Modified: grass/trunk/raster/r.topmodel/infiltration.c
===================================================================
--- grass/trunk/raster/r.topmodel/infiltration.c 2013-10-20 03:57:31 UTC (rev 58072)
+++ grass/trunk/raster/r.topmodel/infiltration.c 2013-10-20 12:19:10 UTC (rev 58073)
@@ -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;
Modified: grass/trunk/raster/r.topmodel/main.c
===================================================================
--- grass/trunk/raster/r.topmodel/main.c 2013-10-20 03:57:31 UTC (rev 58072)
+++ grass/trunk/raster/r.topmodel/main.c 2013-10-20 12:19:10 UTC (rev 58073)
@@ -4,7 +4,7 @@
* TMOD9502.FOR Author: Keith Beven <k.beven at 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 at gmail.com>
* Hydro Laboratory, Kyungpook National University
* South Korea
@@ -14,47 +14,30 @@
*
*/
+#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 Option *topidxclass;
+ struct Option *topidx;
+ struct Option *ntopidxclasses;
+ struct Option *outtopidxstats;
+ } params;
- struct
- {
- struct Flag *input;
- } flag;
-
/* Initialize GRASS and parse command line */
G_gisinit(argv[0]);
@@ -65,154 +48,115 @@
_("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");
+ params.params = G_define_standard_option(G_OPT_F_INPUT);
+ params.params->key = "parameters";
+ params.params->description = _("Name of TOPMODEL parameters file");
- 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");
+ params.topidxstats = G_define_standard_option(G_OPT_F_INPUT);
+ params.topidxstats->key = "topidxstats";
+ params.topidxstats->description =
+ _("Name of topographic index statistics file");
- 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");
+ params.input = G_define_standard_option(G_OPT_F_INPUT);
+ params.input->description =
+ _("Name of rainfall and potential evapotranspiration data file");
- 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");
+ params.output = G_define_standard_option(G_OPT_F_OUTPUT);
+ params.output->description = _("Name for output file");
- 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");
+ 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;
- 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");
+ 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;
- 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 =
- _("Name of rainfall and potential evapotranspiration data file");
+ 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;
- param.output = G_define_standard_option(G_OPT_F_OUTPUT);
- param.output->description = _("Name for output file");
+ 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");
- 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");
+ 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");
- /* Flag definitions */
- flag.input = G_define_flag();
- flag.input->key = 'i';
- flag.input->description = _("Input data given for (o/i)");
+ 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;
+ /* 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;
- misc.nidxclass = atoi(param.nidxclass->answer);
+ if (!params.timestep->answer)
+ params.timestep->answer = "0";
+ misc.timestep = atoi(params.timestep->answer);
- if (!param.timestep->answer)
- param.timestep->answer = "0";
- if (!param.idxclass->answer)
- param.idxclass->answer = "0";
+ if (!params.topidxclass->answer)
+ params.topidxclass->answer = "0";
+ misc.topidxclass = atoi(params.topidxclass->answer);
- misc.timestep = atoi(param.timestep->answer);
- misc.idxclass = atoi(param.idxclass->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;
- flg.input = flag.input->answer;
- flg.overwr = module->overwrite;
+ if (ntopidxclasses <= 0)
+ G_fatal_error(_("ntopidxclasses must be positive."));
- 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();
+ 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);
}
Modified: grass/trunk/raster/r.topmodel/r.topmodel.html
===================================================================
--- grass/trunk/raster/r.topmodel/r.topmodel.html 2013-10-20 03:57:31 UTC (rev 58072)
+++ grass/trunk/raster/r.topmodel/r.topmodel.html 2013-10-20 12:19:10 UTC (rev 58073)
@@ -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 @@
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 @@
<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 at gmail com">Huidae Cho</a><br>
+<a href="mailto:grass4u at 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>
Modified: grass/trunk/raster/r.topmodel/topmodel.c
===================================================================
--- grass/trunk/raster/r.topmodel/topmodel.c 2013-10-20 03:57:31 UTC (rev 58072)
+++ grass/trunk/raster/r.topmodel/topmodel.c 2013-10-20 12:19:10 UTC (rev 58073)
@@ -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];
-double get_lambda(void)
+ sprintf(input, "input=%s", topidx);
+ sprintf(output, "output=%s", outtopidxstats);
+ sprintf(nsteps, "nsteps=%d", ntopidxclasses);
+
+ 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 @@
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 @@
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 @@
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 @@
}
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 @@
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 @@
_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 @@
_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 @@
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 @@
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 @@
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 @@
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 @@
}
}
-
return;
}
-
-void topmodel(void)
+void run_topmodel(void)
{
initialize();
- implement();
- others();
+ calculate_flows();
+ calculate_others();
-
return;
}
More information about the grass-commit
mailing list