[GRASS-SVN] r60790 - grass/branches/releasebranch_7_0/raster/r.topmodel
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jun 11 00:10:46 PDT 2014
Author: hcho
Date: 2014-06-11 00:10:46 -0700 (Wed, 11 Jun 2014)
New Revision: 60790
Modified:
grass/branches/releasebranch_7_0/raster/r.topmodel/file_io.c
grass/branches/releasebranch_7_0/raster/r.topmodel/global.h
grass/branches/releasebranch_7_0/raster/r.topmodel/infiltration.c
grass/branches/releasebranch_7_0/raster/r.topmodel/main.c
grass/branches/releasebranch_7_0/raster/r.topmodel/topmodel.c
Log:
r.topmodel: backport r60788 and r60789 from trunk
Modified: grass/branches/releasebranch_7_0/raster/r.topmodel/file_io.c
===================================================================
--- grass/branches/releasebranch_7_0/raster/r.topmodel/file_io.c 2014-06-11 07:02:46 UTC (rev 60789)
+++ grass/branches/releasebranch_7_0/raster/r.topmodel/file_io.c 2014-06-11 07:10:46 UTC (rev 60790)
@@ -97,7 +97,7 @@
}
if (params.qs0 == 0.0) {
fclose(fp);
- G_fatal_error("parameters.qs0 cannot be 0.0");
+ G_fatal_error(_("%s cannot be 0"), "parameters.qs0");
exit(EXIT_FAILURE);
}
for (; !feof(fp);) {
@@ -233,31 +233,31 @@
ltime->tm_year, ltime->tm_mon, ltime->tm_mday,
ltime->tm_hour, ltime->tm_min, ltime->tm_sec);
fprintf(fp, "#\n");
- fprintf(fp, "# Qt_peak [m^3/timestep]: Total peak flow\n");
+ fprintf(fp, "# Qt_peak [m^3/timestep]: Peak total flow\n");
fprintf(fp, "# tt_peak [timestep]: Peak time for total flow\n");
fprintf(fp, "# Qt_mean [m^3/timestep]: Mean total flow\n");
- fprintf(fp, "# lnTe [ln(m^2/timestep)]: Areal average of ln(T0)\n");
+ fprintf(fp, "# lnTe [ln(m^2/timestep)]: ln of the average transmissivity at the soil surface\n");
fprintf(fp, "# vch [m/timestep]: Main channel routing velocity\n");
fprintf(fp, "# vr [m/timestep]: Internal subcatchment routing velocity\n");
- fprintf(fp, "# lambda [ln(m^2)]: Areal average of topographic index\n");
- fprintf(fp, "# qss [m/timestep]: Subsurface flow per unit area at a soil surface\n");
+ fprintf(fp, "# lambda [ln(m^2)]: Average topographic index\n");
+ fprintf(fp, "# qss [m/timestep]: Saturated subsurface flow per unit area\n");
fprintf(fp, "# qs0 [m/timestep]: Initial subsurface flow per unit area\n");
- fprintf(fp, "# ncells: Number of non-NULL cells\n");
+ fprintf(fp, "# ncells: Number of valid cells\n");
fprintf(fp, "# ntopidxclasses: Number of topographic index classes\n");
fprintf(fp, "# dt [h]: Time step\n");
fprintf(fp, "# ntimesteps: Number of time steps\n");
- fprintf(fp, "# nch: Number of routing time steps\n");
- fprintf(fp, "# nreaches: Number of reach time steps (time of concentration)\n");
- fprintf(fp, "# ndelays: Delay time steps between rainfall and flow response\n");
+ fprintf(fp, "# nch: Number of channel segments\n");
+ fprintf(fp, "# delay [timestep]: Routing delay in the main channel\n");
+ fprintf(fp, "# tc [timestep]: Time of concentration\n");
fprintf(fp, "#\n");
- fprintf(fp, "# tch [timestep]: Routing time\n");
- fprintf(fp, "# Ad [m^2]: Difference in contribution area for reach time steps\n");
+ fprintf(fp, "# tch [timestep]: Routing time to the catchment outlet\n");
+ fprintf(fp, "# Ad [m^2]: Difference in the contribution area\n");
fprintf(fp, "# Qt [m^3/timestep]: Total flow\n");
fprintf(fp, "# qt [m/timestep]: Total flow per unit area\n");
- fprintf(fp, "# qo [m/timestep]: Saturation overland flow per unit area\n");
+ fprintf(fp, "# qo [m/timestep]: Saturated overland flow per unit area\n");
fprintf(fp, "# qs [m/timestep]: Subsurface flow per unit area\n");
- fprintf(fp, "# qv [m/timestep]: Vertical flux (or drainage flux)\n");
- fprintf(fp, "# S_mean [m]: Mean saturation deficit in the watershed\n");
+ fprintf(fp, "# qv [m/timestep]: Vertical drainage flux from unsaturated zone\n");
+ fprintf(fp, "# S_mean [m]: Mean saturation deficit\n");
if (params.infex) {
fprintf(fp, "# f [m/timestep]: Infiltration rate\n");
fprintf(fp, "# fex [m/timestep]: Infiltration excess runoff\n");
@@ -266,10 +266,10 @@
if (misc.timestep || misc.topidxclass) {
fprintf(fp, "#\n");
fprintf(fp, "# Srz [m]: Root zone storage deficit\n");
- fprintf(fp, "# Suz [m]: Unsaturated (gravity drainage) zone storage\n");
+ fprintf(fp, "# Suz [m]: Unsaturated zone storage (gravity drainage)\n");
fprintf(fp, "# S [m]: Local saturated zone deficit due to gravity drainage\n");
fprintf(fp, "# Ea [m/timestep]: Actual evapotranspiration\n");
- fprintf(fp, "# ex [m/timestep]: Excess flow from a fully saturated area per unit area\n");
+ fprintf(fp, "# ex [m/timestep]: Excess flow from fully saturated area per unit area\n");
}
fprintf(fp, "\n");
@@ -287,8 +287,8 @@
fprintf(fp, "dt: %10.3e\n", input.dt);
fprintf(fp, "ntimesteps: %10d\n", input.ntimesteps);
fprintf(fp, "nch: %10d\n", params.nch);
- fprintf(fp, "nreaches: %10d\n", misc.nreaches);
- fprintf(fp, "ndelays: %10d\n", misc.ndelays);
+ fprintf(fp, "delay: %10d\n", misc.delay);
+ fprintf(fp, "tc: %10d\n", misc.tc);
fprintf(fp, "\n");
fprintf(fp, "%10s\n", "tch");
@@ -297,7 +297,7 @@
fprintf(fp, "\n");
fprintf(fp, "%10s\n", "Ad");
- for (i = 0; i < misc.nreaches; i++)
+ for (i = 0; i < misc.tc; i++)
fprintf(fp, "%10.3e\n", misc.Ad[i]);
fprintf(fp, "\n");
Modified: grass/branches/releasebranch_7_0/raster/r.topmodel/global.h
===================================================================
--- grass/branches/releasebranch_7_0/raster/r.topmodel/global.h 2014-06-11 07:02:46 UTC (rev 60789)
+++ grass/branches/releasebranch_7_0/raster/r.topmodel/global.h 2014-06-11 07:10:46 UTC (rev 60790)
@@ -89,8 +89,8 @@
int ncells;
/* Number of topographic index classes */
int ntopidxclasses;
- int ndelays;
- int nreaches;
+ int delay;
+ int tc;
double lnTe;
double vch;
double vr;
@@ -103,7 +103,7 @@
int tt_peak;
/* params.nch's */
double *tch;
- /* misc.nreach's */
+ /* misc.tc's */
double *Ad;
/* input.ntimestep's */
double *Qt;
Modified: grass/branches/releasebranch_7_0/raster/r.topmodel/infiltration.c
===================================================================
--- grass/branches/releasebranch_7_0/raster/r.topmodel/infiltration.c 2014-06-11 07:02:46 UTC (rev 60789)
+++ grass/branches/releasebranch_7_0/raster/r.topmodel/infiltration.c 2014-06-11 07:10:46 UTC (rev 60790)
@@ -66,7 +66,7 @@
}
if (i == MAXITER)
G_warning(
- _("Maximum number of iterations exceeded at timestep %d."),
+ _("Maximum number of iterations exceeded at timestep %d"),
timestep);
pt = t - input.dt + (f - cumf) / R;
@@ -105,7 +105,7 @@
break;
}
if (i == MAXITER)
- G_warning(_("Maximum number of iterations exceeded at timestep %d."),
+ G_warning(_("Maximum number of iterations exceeded at timestep %d"),
timestep);
if (f < cumf + R * input.dt) {
Modified: grass/branches/releasebranch_7_0/raster/r.topmodel/main.c
===================================================================
--- grass/branches/releasebranch_7_0/raster/r.topmodel/main.c 2014-06-11 07:02:46 UTC (rev 60789)
+++ grass/branches/releasebranch_7_0/raster/r.topmodel/main.c 2014-06-11 07:10:46 UTC (rev 60790)
@@ -1,18 +1,21 @@
-/*
- * r.topmodel: simulates TOPMODEL based on TMOD9502.FOR.
+
+/****************************************************************************
*
- * TMOD9502.FOR Author: Keith Beven <k.beven at lancaster.ac.uk>
- * http://www.es.lancs.ac.uk/hfdg/topmodel.html
+ * MODULE: r.topmodel
*
- * Copyright (C) 2000, 2010, 2013 by the GRASS Development Team
- * Author: Huidae Cho <grass4u at gmail.com>
- * Hydro Laboratory, Kyungpook National University
- * South Korea
+ * AUTHOR(S): Huidae Cho <grass4u gmail.com>, Hydro Laboratory,
+ * Kyungpook National University
+ * Based on TMOD9502.FOR by Keith Beven <k.beven lancaster.ac.uk>
*
- * This program is free software under the GPL (>=v2)
- * Read the file COPYING coming with GRASS for details.
+ * PURPOSE: Simulates TOPMODEL.
*
- */
+ * COPYRIGHT: (C) 2000-2014 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
#define _MAIN_C_
#include <stdio.h>
@@ -85,7 +88,7 @@
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");
+ _("Name of input topographic index raster map");
params.topidx->description =
_("Must be clipped to the catchment boundary. Used for generating outtopidxstats.");
params.topidx->required = NO;
@@ -143,13 +146,15 @@
outtopidxstats = params.outtopidxstats->answer;
if (ntopidxclasses <= 0)
- G_fatal_error(_("ntopidxclasses must be positive."));
+ G_fatal_error(_("%s must be positive"), "ntopidxclasses");
create_topidxstats(topidx, ntopidxclasses, outtopidxstats);
} else if (params.topidx->answer) {
- G_warning(_("Ignoring topidx because outtopidxstats is not specified."));
+ G_warning(_("Ignoring %s because %s is not specified"),
+ params.topidx->key, params.outtopidxstats->key);
} else if (params.outtopidxstats->answer) {
- G_warning(_("Ignoring outtopidxstats because topidx is not specified."));
+ G_warning(_("Ignoring %s because %s is not specified"),
+ params.outtopidxstats->key, params.topidx->key);
}
if (flags.preprocess->answer)
Modified: grass/branches/releasebranch_7_0/raster/r.topmodel/topmodel.c
===================================================================
--- grass/branches/releasebranch_7_0/raster/r.topmodel/topmodel.c 2014-06-11 07:02:46 UTC (rev 60789)
+++ grass/branches/releasebranch_7_0/raster/r.topmodel/topmodel.c 2014-06-11 07:10:46 UTC (rev 60790)
@@ -1,6 +1,7 @@
#include <math.h>
#include <grass/raster.h>
#include <grass/spawn.h>
+#include <grass/glocale.h>
#include "global.h"
void create_topidxstats(char *topidx, int ntopidxclasses, char *outtopidxstats)
@@ -17,9 +18,10 @@
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");
+ G_fatal_error(_("Unable to run %s"), "r.stats");
}
+/* Calculate the areal average of topographic index */
double calculate_lambda(void)
{
int i;
@@ -33,36 +35,55 @@
return lambda;
}
+/* Initialize the flows */
void initialize(void)
{
int i, j, t;
double A1, A2;
+ /* average topographic index */
misc.lambda = calculate_lambda();
+
+ /* ln of the average transmissivity at the soil surface */
misc.lnTe = params.lnTe + log(input.dt);
+
+ /* main channel routing velocity */
misc.vch = params.vch * input.dt;
+
+ /* internal subcatchment routing velocity */
misc.vr = params.vr * input.dt;
+
+ /* initial subsurface flow per unit area */
misc.qs0 = params.qs0 * input.dt;
+
+ /* saturated subsurface flow per unit area */
misc.qss = exp(misc.lnTe - misc.lambda);
misc.tch = (double *)G_malloc(params.nch * sizeof(double));
+
+ /* routing time in the main channel */
misc.tch[0] = params.d[0] / misc.vch;
for (i = 1; i < params.nch; i++)
+ /* routing time in each internal subcatchment channel */
misc.tch[i] = misc.tch[0] + (params.d[i] - params.d[0]) / misc.vr;
- 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];
+ /* time of concentration */
+ misc.tc = (int)misc.tch[params.nch - 1];
+ if ((double)misc.tc < misc.tch[params.nch - 1])
+ misc.tc++;
- misc.nreaches -= misc.ndelays;
+ /* routing delay in the main channel */
+ misc.delay = (int)misc.tch[0];
- 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]) {
+ /* time of concentration in the subcatchment */
+ misc.tc -= misc.delay;
+
+ /* cumulative ratio of the contribution area for each time step */
+ misc.Ad = (double *)G_malloc(misc.tc * sizeof(double));
+ for (i = 0; i < misc.tc; i++) {
+ t = misc.delay + i + 1;
+ if (t > misc.tch[params.nch - 1])
misc.Ad[i] = 1.0;
- }
else {
for (j = 1; j < params.nch; j++) {
if (t <= misc.tch[j]) {
@@ -76,14 +97,13 @@
}
}
-
+ /* difference in the contribution area for each time step */
A1 = misc.Ad[0];
misc.Ad[0] *= params.A;
- for (i = 1; i < misc.nreaches; i++) {
+ for (i = 1; i < misc.tc; i++) {
A2 = misc.Ad[i];
- misc.Ad[i] = A2 - A1;
+ misc.Ad[i] = (A2 - A1) * params.A;
A1 = A2;
- misc.Ad[i] *= params.A;
}
misc.Srz = (double **)G_malloc(input.ntimesteps * sizeof(double *));
@@ -94,24 +114,28 @@
}
for (i = 0; i < misc.ntopidxclasses; i++) {
+ /* initial root zone storage deficit */
misc.Srz[0][i] = params.Sr0;
+ /* initial unsaturated zone storage */
misc.Suz[0][i] = 0.0;
}
misc.S_mean = (double *)G_malloc(input.ntimesteps * sizeof(double));
+ /* initial mean saturation deficit */
misc.S_mean[0] = -params.m * log(misc.qs0 / misc.qss);
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.ndelays; i++)
- misc.Qt[i] = misc.qs0 * params.A;
-
+ /* initial total flow */
A1 = 0.0;
- for (i = 0; i < misc.nreaches; i++) {
- A1 += misc.Ad[i];
- misc.Qt[misc.ndelays + i] = misc.qs0 * (params.A - A1);
+ for (i = 0; i < input.ntimesteps; i++) {
+ if (i < misc.delay)
+ misc.Qt[i] = misc.qs0 * params.A;
+ else if (i < misc.delay + misc.tc) {
+ A1 += misc.Ad[i - misc.delay];
+ misc.Qt[i] = misc.qs0 * (params.A - A1);
+ } else
+ misc.Qt[i] = 0.0;
}
}
@@ -119,7 +143,7 @@
{
int i, j, k;
double Aatb_r;
- double R;
+ double f;
double qo, qv;
misc.S = (double **)G_malloc(input.ntimesteps * sizeof(double *));
@@ -152,15 +176,19 @@
misc.qs[i] = 0.0;
if (params.infex) {
+ /* infiltration */
misc.f[i] = input.dt *
calculate_infiltration(i + 1, input.R[i] / input.dt);
+ /* infiltration excess runoff */
misc.fex[i] = input.R[i] - misc.f[i];
- R = misc.f[i];
+ f = misc.f[i];
}
else {
+ /* no infiltration excess runoff */
misc.f[i] = 0.0;
misc.fex[i] = 0.0;
- R = input.R[i];
+ /* 100% of rainfall infiltrates */
+ f = input.R[i];
}
if (i) {
@@ -170,31 +198,40 @@
}
}
+ /* subsurface flow */
misc.qs[i] = misc.qss * exp(-misc.S_mean[i] / params.m);
for (j = 0; j < misc.ntopidxclasses; j++) {
+ /* average area of a topographic index class */
Aatb_r = (topidxstats.Aatb_r[j] +
(j < misc.ntopidxclasses - 1 ? topidxstats.Aatb_r[j + 1]
: 0.0)) / 2.0;
+ /* saturation deficit */
misc.S[i][j] = misc.S_mean[i] +
params.m * (misc.lambda - topidxstats.atb[j]);
if (misc.S[i][j] < 0.0)
+ /* fully saturated */
misc.S[i][j] = 0.0;
- misc.Srz[i][j] -= R;
-
+ /* root zone storage deficit */
+ misc.Srz[i][j] -= f;
if (misc.Srz[i][j] < 0.0) {
+ /* full storage */
+ /* unsaturated zone storage */
misc.Suz[i][j] -= misc.Srz[i][j];
misc.Srz[i][j] = 0.0;
}
+ /* if there is enough unsaturated zone storage */
misc.ex[i][j] = 0.0;
if (misc.Suz[i][j] > misc.S[i][j]) {
+ /* saturation excess */
misc.ex[i][j] = misc.Suz[i][j] - misc.S[i][j];
misc.Suz[i][j] = misc.S[i][j];
}
+ /* drainage from unsaturated zone */
qv = 0.0;
if (misc.S[i][j] > 0.0) {
qv = (params.td > 0.0 ?
@@ -212,6 +249,7 @@
misc.qv[i][j] = qv;
misc.qv[i][misc.ntopidxclasses] += misc.qv[i][j];
+ /* evapotranspiration from root zone storage deficit */
misc.Ea[i][j] = 0.0;
if (input.Ep[i] > 0.0) {
misc.Ea[i][j] = input.Ep[i] *
@@ -221,6 +259,7 @@
}
misc.Srz[i][j] += misc.Ea[i][j];
+ /* overland flow from fully saturated area */
qo = 0.0;
if (j > 0) {
if (misc.ex[i][j] > 0.0)
@@ -234,25 +273,29 @@
misc.qo[i][j] = qo;
misc.qo[i][misc.ntopidxclasses] += misc.qo[i][j];
+ /* total flow */
misc.qt[i][j] = misc.qo[i][j] + misc.qs[i];
}
+ /* aggregate flows over topographic index classes */
misc.qo[i][misc.ntopidxclasses] += misc.fex[i];
- misc.qt[i][misc.ntopidxclasses] = misc.qo[i][misc.ntopidxclasses] + misc.qs[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.ntopidxclasses];
-
+ /* mean saturation deficit */
+ misc.S_mean[i] += misc.qs[i] - misc.qv[i][misc.ntopidxclasses];
if (i + 1 < input.ntimesteps)
misc.S_mean[i + 1] = misc.S_mean[i];
- for (j = 0; j < misc.nreaches; j++) {
- k = i + j + misc.ndelays;
+ /* total flow in m^3/timestep */
+ for (j = 0; j < misc.tc; j++) {
+ k = i + j + misc.delay;
if (k > input.ntimesteps - 1)
break;
misc.Qt[k] += misc.qt[i][misc.ntopidxclasses] * misc.Ad[j];
}
}
+ /* mean total flow */
misc.Qt_mean = 0.0;
for (i = 0; i < input.ntimesteps; i++) {
misc.Qt_mean += misc.Qt[i];
@@ -262,6 +305,9 @@
}
}
misc.Qt_mean /= input.ntimesteps;
+
+ /* time of concentration */
+ misc.tc += misc.delay;
}
void run_topmodel(void)
More information about the grass-commit
mailing list