[GRASS-SVN] r30240 - in
grass/branches/releasebranch_6_3/raster/simwe: r.sim.sediment
r.sim.water simlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Feb 18 15:52:15 EST 2008
Author: neteler
Date: 2008-02-18 15:52:12 -0500 (Mon, 18 Feb 2008)
New Revision: 30240
Modified:
grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/description.html
grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/main.c
grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/description.html
grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/main.c
grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/waterglobs.h
grass/branches/releasebranch_6_3/raster/simwe/simlib/erod.c
grass/branches/releasebranch_6_3/raster/simwe/simlib/hydro.c
grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c
grass/branches/releasebranch_6_3/raster/simwe/simlib/output.c
grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h
Log:
various fixes backported from HEAD
Modified: grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/description.html
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/description.html 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/description.html 2008-02-18 20:52:12 UTC (rev 30240)
@@ -78,7 +78,7 @@
simulation for effective erosion prevention. Water Resources Research,
34(3), 505-516.</p>
<p>
-<a href="http://mpa.itc.it/grasstutor/">
+<a href="http://www.grassbook.org">
Neteler, M. and Mitasova, H., 2004, Open Source GIS: A GRASS GIS Approach, Second Edition, </a>
Kluwer International Series in Engineering and Computer Science, 773, Kluwer Academic Press / Springer,
Boston, Dordrecht, 424 pages.
Modified: grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/main.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/main.c 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/main.c 2008-02-18 20:52:12 UTC (rev 30240)
@@ -1,3 +1,19 @@
+/****************************************************************************
+ *
+ * MODULE: r.simwe.erosion: main program for hydrologic and sediment transport
+ * simulation (SIMWE)
+ *
+ * AUTHOR(S): L. Mitas, H. Mitasova, J. Hofierka
+ * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ *
+ * COPYRIGHT: (C) 2002 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.
+ *
+ *****************************************************************************/
+
/*-
* r.simwe.erosion: main program for hydrologic and sediment transport
* simulation (SIMWE)
Modified: grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/description.html
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/description.html 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/description.html 2008-02-18 20:52:12 UTC (rev 30240)
@@ -9,9 +9,10 @@
resolutions (Mitas and Mitasova 1998). The key inputs of the model include
elevation (<i>elevin</i> raster map), flow gradient vector given by
first-order partial derivatives of elevation field (<i>dxin</i> and <i>dyin</i> raster maps), rainfall
-excess rate (<i>rain</i> raster map) and a surface roughness coefficient
-given by Manning's n (<i>manin</i> raster map). Partial
-derivatives raster maps can be computed along with the interpolation of a DEM using
+excess rate (<i>rain</i> raster map or <i>rain_val</i> single value)
+and a surface roughness coefficient given by Manning's n
+(<i>manin</i> raster map or <i>manin_val</i> single value). Partial
+derivatives raster maps can be computed along with interpolation of a DEM using
the -d option in <a href="v.surf.rst.html">
v.surf.rst</a> module. If elevation raster is already provided, partial derivatives
can be computed using <a href="r.slope.aspect.html">r.slope.aspect</a> module.
@@ -35,7 +36,8 @@
For saturated soil and steady-state water flow it can be estimated using
saturated hydraulic conductivity rates based on field measurements or using
reference values which can be found in literature.
-Optionally, user can provide a runoff infiltration rate map <i>infil</i> in [mm/hr]
+Optionally, user can provide an overland flow infiltration rate map
+<i>infil</i> or a single value <i>infil_val</i> in [mm/hr]
that control the rate of infiltration for the already flowing water, effectively
reducing the flow depth and discharge.
Overland flow can be further controled by permeable check dams or similar type of structures,
@@ -50,16 +52,18 @@
and err is its RMSE). The output vector points map <i>outwalk</i> can be used to analyze and visualize
spatial distribution of walkers at different simulation times (note that
the resulting water depth is based on the density of these walkers). Number
-of the output walkers is controled by the <i>density</i> parameter, which says
-how many walkers used in simalution should be used in the output
+of the output walkers is controled by the <i>density</i> parameter, which controls
+how many walkers used in simulation should be written into the output.
<!--(<font color="#ff0000"> toto treba upresnit/zmenit, lebo nwalk ide prec</font>). -->
Duration of simulation is controled by the <i>niter</i> parameter. The default value
-is 1000 iterations, to reach the steady-state may require, depending on the time step,
-complexity of terrain and land cover and size of the area, several thousand iterations.
-Output files can be saved during simulation using <i>outiter</i> parameter
-defining the time step for writing output files. This option requires
-the time series flag <i>-t</i>. Files are saved with suffix containing
+is 10 minutes, reaching the steady-state may require much longer time,
+depending on the time step, complexity of terrain, land cover and size of the area.
+Output water depth and discharge maps can be saved during simulation using
+the time series flag <i>-t</i> and <i>outiter</i> parameter
+defining the time step in minutes for writing output files.
+Files are saved with suffix containing
iteration number (e.g. name.500, name.1000, etc.).<br>
+<P>
Overland flow is routed based on partial derivatives of elevation
field or other landscape features influencing water flow. Simulation
equations include a diffusion term (<i>diffc</i> parameter) which enables
@@ -145,22 +149,24 @@
Path sampling method for modeling overland water flow, sediment transport
and short term terrain evolution in Open Source GIS.
In: C.T. Miller, M.W. Farthing, V.G. Gray, G.F. Pinder eds.,
-Computational Methods in Water Resources, Elsevier.
+Proceedings of the XVth International Conference on Computational Methods in Water
+Resources (CMWR XV), June 13-17 2004, Chapel Hill, NC, USA, Elsevier, pp. 1479-1490.
<P>
<a href="http://skagit.meas.ncsu.edu/~helena/gmslab/gisc00/duality.html">
- Mitasova H, Mitas, L., 2000, Modeling spatial processes in multiscale framework: exploring duality between particles and fields, plenary talk at GIScience2000 conference, Savannah, GA. </a>
+Mitasova H, Mitas, L., 2000, Modeling spatial processes in multiscale framework:
+exploring duality between particles and fields, </a>
+plenary talk at GIScience2000 conference, Savannah, GA.
<P>
Mitas, L., and Mitasova, H., 1998, Distributed soil erosion simulation
-for effective erosion prevention. Water Resources Research, 34(3), 505-516.<p></p>
+for effective erosion prevention. Water Resources Research, 34(3), 505-516.
<P>
- Mitasova, H., Mitas, L., 2001, Multiscale soil erosion simulations for land use management,
+<a href="http://skagit.meas.ncsu.edu/~helena/gmslab/papers/LLEmiterev1.pdf">
+ Mitasova, H., Mitas, L., 2001, Multiscale soil erosion simulations for land use management, </a>
In: Landscape erosion and landscape evolution modeling, Harmon R. and Doe W. eds.,
Kluwer Academic/Plenum Publishers, pp. 321-347.
<p>
-<a href="http://mpa.itc.it/grasstutor/">
-Neteler, M. and Mitasova, H., 2004, Open Source GIS: A GRASS GIS Approach, Second Edition, </a>
-Kluwer International Series in Engineering and Computer Science, 773, Kluwer Academic Press / Springer,
-Boston, Dordrecht, 424 pages.
+<a href="http://www.grassbook.org">
+Neteler, M. and Mitasova, H., 2008, Open Source GIS: A GRASS GIS Approach. Third Edition.</a>
+The International Series in Engineering and Computer Science: Volume 773. Springer New York Inc, p. 406.
<P>
-Last changed: Date: 2003/11/01 15:55:10 $<p></p>
-</body></html>
+Last changed: Date: 2008/02/16 15:55:10 $<p></p>
Modified: grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/main.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/main.c 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/main.c 2008-02-18 20:52:12 UTC (rev 30240)
@@ -1,5 +1,21 @@
+/****************************************************************************
+ *
+ * MODULE: r.sim.water: main program for hydrologic and sediment transport
+ * simulation (SIMWE)
+ *
+ * AUTHOR(S): L. Mitas, H. Mitasova, J. Hofierka
+ * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ *
+ * COPYRIGHT: (C) 2002 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.
+ *
+ *****************************************************************************/
+
/*-
- * r.simwe.hydro: main program for hydrologic and sediment transport
+ * r.sim.water: main program for hydrologic and sediment transport
* simulation (SIMWE)
*
* Original program (2002) and various modifications:
@@ -27,18 +43,27 @@
*
* Notes on modifications:
* v. 1.0 May 2002
- *
+ * modified by Y. Chemin in February 2008 (reporting, optional inputs)
*/
-#define NWALK "2000000"
+/********************************/
+/* DEFINE GLOB VAR */
+/********************************/
+/* #define NWALK "1000000" */
#define DIFFC "0.8"
-#define HMAX "0.4"
+#define HMAX "0.3"
#define HALPHA "4.0"
#define HBETA "0.5"
-#define NITER "1200"
-#define ITEROUT "300"
+#define NITER "10"
+#define ITEROUT "2"
#define DENSITY "200"
+#define RAINVAL "50"
+#define MANINVAL "0.1"
+#define INFILVAL "0.0"
+/********************************/
+/* INCLUDES */
+/********************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
@@ -54,6 +79,9 @@
#include <grass/site.h>
#include <grass/glocale.h>
+/********************************/
+/* Specific stuff */
+/********************************/
#define MAIN
#include <grass/waterglobs.h>
@@ -67,352 +95,434 @@
char msg[1024];
+/****************************************/
+/* MAIN */
+/****************************************/
int main ( int argc, char *argv[])
{
- int i,ii,j,l;
- int ret_val;
- double x_orig, y_orig;
- static int rand1 = 12345;
- static int rand2 = 67891;
+ int i,ii,j,l;
+ int ret_val;
+ double x_orig, y_orig;
+ static int rand1 = 12345;
+ static int rand2 = 67891;
+
+ G_gisinit (argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("raster, flow, hydrology");
+ module->description =
+ _("Overland flow hydrologic simulation using "
+ "path sampling method (SIMWE)");
+
+ if (G_get_set_window (&cellhd) == -1)
+ exit (EXIT_FAILURE);
- G_gisinit (argv[0]);
+ conv = G_database_units_to_meters_factor();
+
+ mixx = conv * cellhd.west;
+ maxx = conv * cellhd.east;
+ miyy = conv * cellhd.south;
+ mayy = conv * cellhd.north;
+
+ stepx = cellhd.ew_res * conv;
+ stepy = cellhd.ns_res * conv;
+ /* step = amin1(stepx,stepy);*/
+ step = (stepx + stepy) / 2.;
+ mx = cellhd.cols;
+ my = cellhd.rows;
+ x_orig = cellhd.west * conv;
+ y_orig = cellhd.south * conv;/* do we need this? */
+ xmin = 0.;
+ ymin = 0.;
+ xp0 = xmin + stepx / 2.;
+ yp0 = ymin + stepy / 2.;
+ xmax = xmin + stepx * (float) mx;
+ ymax = ymin + stepy * (float) my;
+ /*
+ bxmi=2093113. * conv;
+ bymi=731331. * conv;
+ bxma=2093461. * conv;
+ byma=731529. * conv;
+ bresx=2. * conv;
+ bresy=2. * conv;
+ maxwab=100000;
+
+ mx2o= (int)((bxma-bxmi)/bresx);
+ my2o= (int)((byma-bymi)/bresy);
+ */
+ /* relative small box coordinates: leave 1 grid layer for overlap */
+ /*
+ bxmi = bxmi - mixx + stepx;
+ bymi = bymi - miyy + stepy;
+ bxma = bxma - mixx - stepx;
+ byma = byma - miyy - stepy;
+ mx2 = mx2o - 2*((int) (stepx / bresx));
+ my2 = my2o - 2*((int) (stepy / bresy));
+ */
+ parm.elevin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.elevin->key = "elevin";
+ parm.elevin->description = _("Name of the elevation raster map [m]");
+ parm.elevin->guisection = _("Input_options");
+
+ parm.dxin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.dxin->key = "dxin";
+ parm.dxin->description = _("Name of the x-derivatives raster map [m/m]");
+ parm.dxin->guisection = _("Input_options");
+
+ parm.dyin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.dyin->key = "dyin";
+ parm.dyin->description = _("Name of the y-derivatives raster map [m/m]");
+ parm.dyin->guisection = _("Input_options");
+
+ parm.rain = G_define_standard_option(G_OPT_R_INPUT);
+ parm.rain->key = "rain";
+ parm.rain->required = NO;
+ parm.rain->description = _("Name of the rainfall excess rate (rain-infilt) raster map [mm/hr]");
+ parm.rain->guisection = _("Input_options");
+
+ parm.rainval = G_define_option();
+ parm.rainval->key = "rain_val";
+ parm.rainval->type = TYPE_DOUBLE;
+ parm.rainval->answer = RAINVAL;
+ parm.rainval->required = NO;
+ parm.rainval->description = _("Rainfall excess rate unique value [mm/hr]");
+ parm.rainval->guisection = _("Input_options");
+
+ parm.infil = G_define_standard_option(G_OPT_R_INPUT);
+ parm.infil->key = "infil";
+ parm.infil->required = NO;
+ parm.infil->description = _("Name of the runoff infiltration rate raster map [mm/hr]");
+ parm.infil->guisection = _("Input_options");
+
+ parm.infilval = G_define_option();
+ parm.infilval->key = "infil_val";
+ parm.infilval->type = TYPE_DOUBLE;
+ parm.infilval->answer = INFILVAL;
+ parm.infilval->required = NO;
+ parm.infilval->description = _("Runoff infiltration rate unique value [mm/hr]");
+ parm.infilval->guisection = _("Input_options");
+
+ parm.manin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.manin->key = "manin";
+ parm.manin->required = NO;
+ parm.manin->description = _("Name of the Mannings n raster map");
+ parm.manin->guisection = _("Input_options");
+
+ parm.maninval = G_define_option();
+ parm.maninval->key = "manin_val";
+ parm.maninval->type = TYPE_DOUBLE;
+ parm.maninval->answer = MANINVAL;
+ parm.maninval->required = NO;
+ parm.maninval->description = _("Mannings n unique value");
+ parm.maninval->guisection = _("Input_options");
+
+ parm.traps = G_define_standard_option(G_OPT_R_INPUT);
+ parm.traps->key = "traps";
+ parm.traps->required = NO;
+ parm.traps->description = _("Name of the flow controls raster map (permeability ratio 0-1)");
+ parm.traps->guisection = _("Input_options");
+
+ /* needs to be updated to GRASS 6 vector format !! */
+ parm.sfile = G_define_standard_option (G_OPT_V_INPUT);
+ parm.sfile->key = "vector";
+ parm.sfile->required = NO;
+ parm.sfile->description = _("Name of the sampling locations vector points map");
+ parm.sfile->guisection = _("Input_options");
+
+ parm.depth = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.depth->key = "depth";
+ parm.depth->required = NO;
+ parm.depth->description = _("Output water depth raster map [m]");
+ parm.depth->guisection = _("Output_options");
+
+ parm.disch = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.disch->key = "disch";
+ parm.disch->required = NO;
+ parm.disch->description = _("Output water discharge raster map [m3/s]");
+ parm.disch->guisection = _("Output_options");
+
+ parm.err = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.err->key = "err";
+ parm.err->required = NO;
+ parm.err->description = _("Output simulation error raster map [m]");
+ parm.err->guisection = _("Output_options");
+
+ parm.outwalk = G_define_standard_option(G_OPT_V_OUTPUT);
+ parm.outwalk->key = "outwalk";
+ parm.outwalk->required = NO;
+ parm.outwalk->description = _("Name of the output walkers vector points map");
+ parm.outwalk->guisection = _("Output_options");
+
+ parm.nwalk = G_define_option();
+ parm.nwalk->key = "nwalk";
+ parm.nwalk->type = TYPE_INTEGER;
+/* parm.nwalk->answer = NWALK; */
+ parm.nwalk->required = NO;
+ parm.nwalk->description = _("Number of walkers, default is twice the no. of cells");
+ parm.nwalk->guisection = _("Parameters");
+
+ parm.niter = G_define_option();
+ parm.niter->key = "niter";
+ parm.niter->type = TYPE_INTEGER;
+ parm.niter->answer = NITER;
+ parm.niter->required = NO;
+ parm.niter->description = _("Time used for iterations [minutes]");
+ parm.niter->guisection = _("Parameters");
+
+ parm.outiter = G_define_option();
+ parm.outiter->key = "outiter";
+ parm.outiter->type = TYPE_INTEGER;
+ parm.outiter->answer = ITEROUT;
+ parm.outiter->required = NO;
+ parm.outiter->description = _("Time interval for creating output maps [minutes]");
+ parm.outiter->guisection = _("Parameters");
+
+ parm.density = G_define_option();
+ parm.density->key = "density";
+ parm.density->type = TYPE_INTEGER;
+ parm.density->answer = DENSITY;
+ parm.density->required = NO;
+ parm.density->description = _("Density of output walkers");
+ parm.density->guisection = _("Parameters");
+
+ parm.diffc = G_define_option();
+ parm.diffc->key = "diffc";
+ parm.diffc->type = TYPE_DOUBLE;
+ parm.diffc->answer = DIFFC;
+ parm.diffc->required = NO;
+ parm.diffc->description = _("Water diffusion constant");
+ parm.diffc->guisection = _("Parameters");
+
+ parm.hmax = G_define_option();
+ parm.hmax->key = "hmax";
+ parm.hmax->type = TYPE_DOUBLE;
+ parm.hmax->answer = HMAX;
+ parm.hmax->required = NO;
+ parm.hmax->description = _("Threshold water depth [m] (diffusion increases after this water depth is reached)");
+ parm.hmax->guisection = _("Parameters");
+
+ parm.halpha = G_define_option();
+ parm.halpha->key = "halpha";
+ parm.halpha->type = TYPE_DOUBLE;
+ parm.halpha->answer = HALPHA;
+ parm.halpha->required = NO;
+ parm.halpha->description = _("Diffusion increase constant");
+ parm.halpha->guisection = _("Parameters");
+
+ parm.hbeta = G_define_option();
+ parm.hbeta->key = "hbeta";
+ parm.hbeta->type = TYPE_DOUBLE;
+ parm.hbeta->answer = HBETA;
+ parm.hbeta->required = NO;
+ parm.hbeta->description = _("Weighting factor for water flow velocity vector");
+ parm.hbeta->guisection = _("Parameters");
+ /*
+ flag.mscale = G_define_flag ();
+ flag.mscale->key = 'm';
+ flag.mscale->description = _("Multiscale simulation");
+ */
+ flag.tserie = G_define_flag ();
+ flag.tserie->key = 't';
+ flag.tserie->description = _("Time-series output");
- module = G_define_module();
- module->keywords = _("raster");
- module->description =
- _("Overland flow hydrologic model based on duality "
- "particle-field concept (SIMWE)");
-
- if (G_get_set_window (&cellhd) == -1)
- exit (EXIT_FAILURE);
- conv = G_database_units_to_meters_factor();
+ if (G_parser (argc, argv))
+ exit (EXIT_FAILURE);
- mixx = conv * cellhd.west;
- maxx = conv * cellhd.east;
- miyy = conv * cellhd.south;
- mayy = conv * cellhd.north;
-
- stepx = cellhd.ew_res * conv;
- stepy = cellhd.ns_res * conv;
-/* step = amin1(stepx,stepy);*/
- step = (stepx + stepy) / 2.;
- mx = cellhd.cols;
- my = cellhd.rows;
- x_orig = cellhd.west * conv;
- y_orig = cellhd.south * conv;/* do we need this? */
- xmin = 0.;
- ymin = 0.;
- xp0 = xmin + stepx / 2.;
- yp0 = ymin + stepy / 2.;
- xmax = xmin + stepx * (float) mx;
- ymax = ymin + stepy * (float) my;
+ /* mscale=flag.mscale->answer;*/
+ ts=flag.tserie->answer;
- bxmi=2093113. * conv;
- bymi=731331. * conv;
- bxma=2093461. * conv;
- byma=731529. * conv;
- bresx=2. * conv;
- bresy=2. * conv;
- maxwab=100000;
+ elevin = parm.elevin->answer;
+ dxin = parm.dxin->answer;
+ dyin = parm.dyin->answer;
+ /* maskmap = parm.maskmap->answer;*/
+ rain = parm.rain->answer;
+ infil = parm.infil->answer;
+ traps = parm.traps->answer;
+ manin = parm.manin->answer;
+ depth = parm.depth->answer;
+ disch = parm.disch->answer;
+ err = parm.err->answer;
+ outwalk = parm.outwalk->answer;
+ sfile = parm.sfile->answer;
+
+/* sscanf(parm.nwalk->answer, "%d", &maxwa); */
+ sscanf(parm.niter->answer, "%d", ×ec);
+ sscanf(parm.outiter->answer, "%d", &iterout);
+ sscanf(parm.density->answer, "%d", &ldemo);
+ sscanf(parm.diffc->answer, "%lf", &frac);
+ sscanf(parm.hmax->answer, "%lf", &hhmax);
+ sscanf(parm.halpha->answer, "%lf", &halpha);
+ sscanf(parm.hbeta->answer, "%lf", &hbeta);
- mx2o= (int)((bxma-bxmi)/bresx);
- my2o= (int)((byma-bymi)/bresy);
+ /* if no rain map input, then:*/
+ if(parm.rain->answer==NULL){
+ /*Check for Rain Unique Value Input*/
+ /* if no rain unique value input*/
+ if(parm.rainval->answer==NULL){
+ /*No rain input so use default*/
+ sscanf(RAINVAL, "%lf", &rain_val);
+ /* if rain unique input exist, load it*/
+ } else {
+ /*Unique value input only*/
+ sscanf(parm.rainval->answer, "%lf", &rain_val);
+ }
+ /* if Rain map exists*/
+ } else {
+ /*Map input, so set rain_val to -999.99*/
+ if(parm.rainval->answer==NULL){
+ rain_val = -999.99;
+ } else {
+ /*both map and unique value exist*/
+ /*Choose the map, discard the unique value*/
+ rain_val = -999.99;
+ }
+ }
+ /* Report the final value of rain_val*/
+ G_message(_("rain_val is set to: %f\n"),rain_val);
-/* relative small box coordinates: leave 1 grid layer for overlap */
+ /* if no Mannings map, then:*/
+ if(parm.manin->answer==NULL){
+ /*Check for Manin Unique Value Input*/
+ /* if no Mannings unique value input*/
+ if(parm.maninval->answer==NULL){
+ /*No Mannings input so use default*/
+ sscanf(MANINVAL, "%lf", &manin_val);
+ /* if mannings unique input value exists, load it*/
+ } else {
+ /*Unique value input only*/
+ sscanf(parm.maninval->answer, "%lf", &manin_val);
+ }
+ /* if Mannings map exists*/
+ } else {
+ /* Map input, set manin_val to -999.99*/
+ if(parm.maninval->answer==NULL){
+ manin_val = -999.99;
+ } else {
+ /*both map and unique value exist*/
+ /*Choose map, discard the unique value*/
+ manin_val = -999.99;
+ }
+ }
+ /* Report the final value of manin_val*/
+ G_message(_("manin_val is set to: %f\n"),manin_val);
- bxmi = bxmi - mixx + stepx;
- bymi = bymi - miyy + stepy;
- bxma = bxma - mixx - stepx;
- byma = byma - miyy - stepy;
- mx2 = mx2o - 2*((int) (stepx / bresx));
- my2 = my2o - 2*((int) (stepy / bresy));
+ /* if no infiltration map, then:*/
+ if(parm.infil->answer==NULL){
+ /*Check for Infil Unique Value Input*/
+ /*if no infiltration unique value input*/
+ if(parm.infilval->answer==NULL){
+ /*No infiltration unique value so use default*/
+ sscanf(INFILVAL, "%lf", &infil_val);
+ /* if infiltration unique value exists, load it*/
+ } else {
+ /*unique value input only*/
+ sscanf(parm.infilval->answer, "%lf", &infil_val);
+ }
+ /* if infiltration map exists */
+ } else {
+ /* Map input, set infil_val to -999.99 */
+ if(parm.infilval->answer==NULL){
+ infil_val = -999.99;
+ } else {
+ /*both map and unique value exist*/
+ /*Choose map, discard the unique value*/
+ infil_val = -999.99;
+ }
+ }
+ /* Report the final value of infil_val*/
+ G_message(_("infil_val is set to: %f\n"),infil_val);
- parm.elevin = G_define_option();
- parm.elevin->key = "elevin";
- parm.elevin->type = TYPE_STRING;
- parm.elevin->required = YES;
- parm.elevin->gisprompt = "old,cell,raster";
- parm.elevin->description = _("Name of the elevation raster map, same units as N,E: [m] or [ft] ");
- parm.elevin->guisection = _("Input_options");
+ /* Recompute timesec from user input in minutes
+ * to real timesec in seconds */
+ timesec = timesec * 60.0;
+ iterout = iterout * 60.0;
+ if ((timesec/iterout) > 100.0)
+ G_message(_("More than 100 files are going to be created !!!!!"));
- parm.dxin = G_define_option();
- parm.dxin->key = "dxin";
- parm.dxin->type = TYPE_STRING;
- parm.dxin->required = YES;
- parm.dxin->gisprompt = "old,cell,raster";
- parm.dxin->description = _("Name of the x-derivatives raster map [m/m] or [ft/ft]");
- parm.dxin->guisection = _("Input_options");
+ /* compute how big the raster is and set this to appr 2 walkers per cell */
+ if(parm.nwalk->answer==NULL){
+ maxwa = mx*my*2;
+ rwalk = (double) (mx*my*2.);
+ G_message(_("default nwalk=%d, rwalk=%f"), maxwa, rwalk);
+ } else {
+ sscanf(parm.nwalk->answer, "%d", &maxwa);
+ rwalk = (double) maxwa;
+ }
- parm.dyin = G_define_option();
- parm.dyin->key = "dyin";
- parm.dyin->type = TYPE_STRING;
- parm.dyin->required = YES;
- parm.dyin->gisprompt = "old,cell,raster";
- parm.dyin->description = _("Name of the y-derivatives raster map [m/m] or [ft/ft]");
- parm.dyin->guisection = _("Input_options");
+/* rwalk = (double) maxwa; */
- parm.rain = G_define_option();
- parm.rain->key = "rain";
- parm.rain->type = TYPE_STRING;
- parm.rain->required = YES;
- parm.rain->gisprompt = "old,cell,raster";
- parm.rain->description = _("Name of the rainfall excess rate (rain-infilt) raster map [mm/hr]");
- parm.rain->guisection = _("Input_options");
+ if (conv != 1.0)
+ G_message(_("Using metric conversion factor %f, step=%f"), conv, step);
- parm.infil = G_define_option();
- parm.infil->key = "infil";
- parm.infil->type = TYPE_STRING;
- parm.infil->required = YES;
- parm.infil->gisprompt = "old,cell,raster";
- parm.infil->description = _("Name of the runoff infiltration rate raster map [mm/hr]");
- parm.infil->guisection = _("Input_options");
+ /*
+ * G_set_embedded_null_value_mode(1);
+ */
- parm.traps = G_define_option();
- parm.traps->key = "traps";
- parm.traps->type = TYPE_STRING;
- parm.traps->required = NO;
- parm.traps->gisprompt = "old,cell,raster";
- parm.traps->description = _("Name of the flow controls raster map (permeability ratio 0-1)");
- parm.traps->guisection = _("Input_options");
+ if ((depth == NULL) && (disch == NULL) && (err == NULL) && (outwalk == NULL))
+ G_warning(_("You are not outputting any raster or vector points maps"));
+ ret_val = input_data();
+ if (ret_val != 1)
+ G_fatal_error(_("Input failed"));
- parm.manin = G_define_option();
- parm.manin->key = "manin";
- parm.manin->type = TYPE_STRING;
- parm.manin->required = YES;
- parm.manin->gisprompt = "old,cell,raster";
- parm.manin->description = _("Name of the Mannings n raster map");
- parm.manin->guisection = _("Input_options");
+ /* memory allocation for output grids */
-/* needs to be updated to GRASS 6 vector format !! */
- parm.sfile = G_define_option ();
- parm.sfile->key = "vector";
- parm.sfile->type = TYPE_STRING;
- parm.sfile->required = NO;
- parm.sfile->gisprompt = "old,vector,vector";
- parm.sfile->description = _("Name of the vector points map with x,y locations");
- parm.sfile->guisection = _("Input_options");
+ gama = (double **)G_malloc (sizeof(double *)*(my));
+ for(l=0;l<my;l++)
+ {
+ gama[l] = (double*)G_malloc (sizeof(double)*(mx));
+ }
+ for (j = 0; j < my; j++)
+ {
+ for (i = 0; i < mx; i++)
+ gama[j][i] = 0.;
+ }
- parm.depth = G_define_option();
- parm.depth->key = "depth";
- parm.depth->type = TYPE_STRING;
- parm.depth->required = NO;
- parm.depth->gisprompt = "new,cell,raster";
- parm.depth->description = _("Output water depth raster map");
- parm.depth->guisection = _("Output_options");
+ if (err != NULL)
+ {
+ gammas = (double **)G_malloc (sizeof(double *)*(my));
+ for(l=0;l<my;l++)
+ {
+ gammas[l] = (double*)G_malloc (sizeof(double)*(mx));
+ }
+ for (j = 0; j < my; j++)
+ {
+ for (i = 0; i < mx; i++)
+ gammas[j][i] = 0.;
+ }
+ }
- parm.disch = G_define_option();
- parm.disch->key = "disch";
- parm.disch->type = TYPE_STRING;
- parm.disch->required = NO;
- parm.disch->gisprompt = "new,cell,raster";
- parm.disch->description = _("Output water discharge raster map");
- parm.disch->guisection = _("Output_options");
+ dif = (float **)G_malloc (sizeof(float *)*(my));
+ for(l=0;l<my;l++)
+ {
+ dif[l] = (float*)G_malloc (sizeof(float)*(mx));
+ }
+ for (j = 0; j < my; j++)
+ {
+ for (i = 0; i < mx; i++)
+ dif[j][i] = 0.;
+ }
+
+ /* if (maskmap != NULL)
+ bitmask = BM_create (cols, rows);
+ IL_create_bitmask (¶ms, bitmask);
+ */
+ seeds(rand1,rand2);
+ grad_check();
+ main_loop();
- parm.err = G_define_option();
- parm.err->key = "err";
- parm.err->type = TYPE_STRING;
- parm.err->required = NO;
- parm.err->gisprompt = "new,cell,raster";
- parm.err->description = _("Output simulation error raster map");
- parm.err->guisection = _("Output_options");
+ if (ts == 0) {
+ ii=output_data (0,1.);
+ if (ii != 1)
+ G_fatal_error(_("Cannot write raster maps"));
+ }
- parm.outwalk = G_define_option ();
- parm.outwalk->key = "outwalk";
- parm.outwalk->type = TYPE_STRING;
- parm.outwalk->required = NO;
- parm.outwalk->gisprompt = "new,vector,vector";
- parm.outwalk->description = _("Name of the output walkers vector points map");
- parm.outwalk->guisection = _("Output_options");
+ if (fdwalkers != NULL)
+ fclose (fdwalkers);
- parm.nwalk = G_define_option();
- parm.nwalk->key = "nwalk";
- parm.nwalk->type = TYPE_INTEGER;
- parm.nwalk->answer = NWALK;
- parm.nwalk->required = NO;
- parm.nwalk->description = _("Number of walkers");
- parm.nwalk->guisection = _("Parameters");
-
- parm.niter = G_define_option();
- parm.niter->key = "niter";
- parm.niter->type = TYPE_INTEGER;
- parm.niter->answer = NITER;
- parm.niter->required = NO;
- parm.niter->description = _("Number of time iterations (sec.)");
- parm.niter->guisection = _("Parameters");
-
- parm.outiter = G_define_option();
- parm.outiter->key = "outiter";
- parm.outiter->type = TYPE_INTEGER;
- parm.outiter->answer = ITEROUT;
- parm.outiter->required = NO;
- parm.outiter->description = _("Time step for saving output maps (sec.)");
- parm.outiter->guisection = _("Parameters");
-
- parm.density = G_define_option();
- parm.density->key = "density";
- parm.density->type = TYPE_INTEGER;
- parm.density->answer = DENSITY;
- parm.density->required = NO;
- parm.density->description = _("Density of output walkers");
- parm.density->guisection = _("Parameters");
-
- parm.diffc = G_define_option();
- parm.diffc->key = "diffc";
- parm.diffc->type = TYPE_DOUBLE;
- parm.diffc->answer = DIFFC;
- parm.diffc->required = NO;
- parm.diffc->description = _("Water diffusion constant");
- parm.diffc->guisection = _("Parameters");
-
- parm.hmax = G_define_option();
- parm.hmax->key = "hmax";
- parm.hmax->type = TYPE_DOUBLE;
- parm.hmax->answer = HMAX;
- parm.hmax->required = NO;
- parm.hmax->description = _("Threshold water depth (diffusion increases after this water depth is reached)");
- parm.hmax->guisection = _("Parameters");
-
- parm.halpha = G_define_option();
- parm.halpha->key = "halpha";
- parm.halpha->type = TYPE_DOUBLE;
- parm.halpha->answer = HALPHA;
- parm.halpha->required = NO;
- parm.halpha->description = _("Diffusion increase constant");
- parm.halpha->guisection = _("Parameters");
-
- parm.hbeta = G_define_option();
- parm.hbeta->key = "hbeta";
- parm.hbeta->type = TYPE_DOUBLE;
- parm.hbeta->answer = HBETA;
- parm.hbeta->required = NO;
- parm.hbeta->description = _("Weighting factor for water flow velocity vector");
- parm.hbeta->guisection = _("Parameters");
-
- flag.mscale = G_define_flag ();
- flag.mscale->key = 'm';
- flag.mscale->description = _("Multiscale simulation");
-
- flag.tserie = G_define_flag ();
- flag.tserie->key = 't';
- flag.tserie->description = _("Time-series (dynamic) output");
-
-
- if (G_parser (argc, argv))
- exit (EXIT_FAILURE);
-
- mscale=flag.mscale->answer;
- ts=flag.tserie->answer;
-
- elevin = parm.elevin->answer;
- dxin = parm.dxin->answer;
- dyin = parm.dyin->answer;
-/* maskmap = parm.maskmap->answer;*/
- rain = parm.rain->answer;
- infil = parm.infil->answer;
- traps = parm.traps->answer;
- manin = parm.manin->answer;
- depth = parm.depth->answer;
- disch = parm.disch->answer;
- err = parm.err->answer;
- outwalk = parm.outwalk->answer;
- sfile = parm.sfile->answer;
-
- sscanf(parm.nwalk->answer, "%d", &maxwa);
- sscanf(parm.niter->answer, "%d", ×ec);
- sscanf(parm.outiter->answer, "%d", &iterout);
- sscanf(parm.density->answer, "%d", &ldemo);
- sscanf(parm.diffc->answer, "%lf", &frac);
- sscanf(parm.hmax->answer, "%lf", &hhmax);
- sscanf(parm.halpha->answer, "%lf", &halpha);
- sscanf(parm.hbeta->answer, "%lf", &hbeta);
-
- /* compute how big the raster is and set this to appr 2 walkers per cell */
- rwalk = (double) maxwa;
-
- if (conv != 1.0)
- G_message(_("Using metric conversion factor %f, step=%f"), conv, step);
-
-
- /*
- * G_set_embedded_null_value_mode(1);
- */
-
- if ((depth == NULL) && (disch == NULL) && (err == NULL) && (outwalk == NULL))
- G_warning(_("You are not outputting any raster or vector points maps"));
-
- ret_val = input_data();
- if (ret_val != 1)
- G_fatal_error(_("Input failed"));
-
-/* memory allocation for output grids */
-
- gama = (double **)G_malloc (sizeof(double *)*(my));
- for(l=0;l<my;l++)
- {
- gama[l] = (double*)G_malloc (sizeof(double)*(mx));
- }
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- gama[j][i] = 0.;
- }
-
- if (err != NULL)
- {
- gammas = (double **)G_malloc (sizeof(double *)*(my));
- for(l=0;l<my;l++)
- {
- gammas[l] = (double*)G_malloc (sizeof(double)*(mx));
- }
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- gammas[j][i] = 0.;
- }
- }
-
- dif = (float **)G_malloc (sizeof(float *)*(my));
- for(l=0;l<my;l++)
- {
- dif[l] = (float*)G_malloc (sizeof(float)*(mx));
- }
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- dif[j][i] = 0.;
- }
-
-
-
-
-/* if (maskmap != NULL)
- bitmask = BM_create (cols, rows);
-
- IL_create_bitmask (¶ms, bitmask);
-*/
-
-
- seeds(rand1,rand2);
-
- grad_check();
-
-/* fprintf (stdout, "\n Percent complete: ");*/
-
- main_loop();
-
- if (ts == 0) {
- ii=output_data (0,1.);
- if (ii != 1)
- G_fatal_error(_("Cannot write raster maps"));
- }
-
- if (fdwalkers != NULL)
- fclose (fdwalkers);
-
- if (sfile != NULL)
- fclose(fw);
-
- exit (EXIT_SUCCESS);
+ if (sfile != NULL)
+ fclose(fw);
+ /* Exit with Success */
+ exit (EXIT_SUCCESS);
}
Modified: grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/waterglobs.h
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/waterglobs.h 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/waterglobs.h 2008-02-18 20:52:12 UTC (rev 30240)
@@ -20,7 +20,8 @@
char *err;
char *outwalk;
char *mapset;
-char *mscale;
+
+/* Flag*/
char *tserie;
char *wdepth;
@@ -33,16 +34,21 @@
char *flux;
char *erdep;
+/*New Input val*/
+char *rainval;
+char *maninval;
+char *infilval;
+
struct
{
struct Option *elevin,*dxin,*dyin,*rain,*infil,*traps,*manin,*sfile,*depth,*disch,*err,
*outwalk,*nwalk,*niter,*outiter,*density,*diffc,*hmax,*halpha,*hbeta,*wdepth,
-*detin,*tranin,*tauin,*tc,*et,*conc,*flux,*erdep;
+*detin,*tranin,*tauin,*tc,*et,*conc,*flux,*erdep,*rainval,*maninval,*infilval;
} parm;
struct
{
- struct Flag *mscale, *tserie;
+ struct Flag *tserie;
} flag;
@@ -82,7 +88,7 @@
int mx, my;
int mx2, my2;
-double bxmi,bymi,bxma,byma,bresx,bresy;
+/*double bxmi,bymi,bxma,byma,bresx,bresy;*/
int maxwab;
double step,conv;
@@ -112,6 +118,11 @@
int miter,nwalka,lwwfin;
double timec;
int ts, timesec;
+
+double rain_val;
+double manin_val;
+double infil_val;
+
#else
extern FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps, *fdmanin, *fddepth, *fddisch, *fderr, *fdoutwalk,*fdwalkers;
extern FILE *fdwdepth,*fddetin,*fdtranin,*fdtauin, *fdtc, *fdet, *fdconc, *fdflux,*fderdep;
@@ -130,7 +141,6 @@
extern char *err;
extern char *outwalk;
extern char *mapset;
-extern char *mscale;
extern char *tserie;
extern char *wdepth;
@@ -143,17 +153,20 @@
extern char *flux;
extern char *erdep;
+extern char *rainval;
+extern char *maninval;
+extern char *infilval;
extern struct
{
struct Option *elevin,*dxin,*dyin,*rain,*infil,*traps,*manin,*sfile,*depth,*disch,*err,
*outwalk,*nwalk,*niter,*outiter,*density,*diffc,*hmax,*halpha,*hbeta,*wdepth,
-*detin,*tranin,*tauin,*tc,*et,*conc,*flux,*erdep;
+*detin,*tranin,*tauin,*tc,*et,*conc,*flux,*erdep,*rainval,*maninval,*infilval;
} parm;
extern struct
{
- struct Flag *mscale, *tserie;
+ struct Flag *tserie;
} flag;
@@ -194,7 +207,7 @@
extern int mx, my;
extern int mx2, my2;
-extern double bxmi,bymi,bxma,byma,bresx,bresy;
+/*extern double bxmi,bymi,bxma,byma,bresx,bresy;*/
extern int maxwab;
extern double step,conv;
@@ -224,4 +237,8 @@
extern int miter,nwalka,lwwfin;
extern double timec;
extern int ts, timesec;
+
+extern double rain_val;
+extern double manin_val;
+extern double infil_val;
#endif
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/erod.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/erod.c 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/erod.c 2008-02-18 20:52:12 UTC (rev 30240)
@@ -1,3 +1,17 @@
+/****************************************************************************
+ *
+ * MODULE: simwe library
+ * AUTHOR(S): Helena Mitasova, Jaro Hofierka, Lubos Mitas:
+ * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ *
+ * COPYRIGHT: (C) 2002 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.
+ *
+ *****************************************************************************/
+
/* erod.c (simlib), 20.nov.2002, JH */
#include <stdio.h>
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/hydro.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/hydro.c 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/hydro.c 2008-02-18 20:52:12 UTC (rev 30240)
@@ -1,3 +1,17 @@
+/****************************************************************************
+ *
+ * MODULE: simwe library
+ * AUTHOR(S): Helena Mitasova, Jaro Hofierka, Lubos Mitas:
+ * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ *
+ * COPYRIGHT: (C) 2002 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.
+ *
+ *****************************************************************************/
+
/* hydro.c (simlib), 20.nov.2002, JH */
#include <stdio.h>
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c 2008-02-18 20:52:12 UTC (rev 30240)
@@ -25,664 +25,720 @@
int input_data (void)
{
- FCELL *cell1, *cell4b,*cell5;
- FCELL *cell9, *cell10,*cell11;
- DCELL *cell2, *cell3, *cell4, *cell4a,*cell12;
- int fd1, fd2, fd3, fd4, fd4a, fd4b, fd5, row, row_rev;
- int fd9, fd10, fd11, fd12;
- int l,j;
- int nn,cc,ii,dd;
- double unitconv=0.0000002; /* mm/hr to m/s */
- char *mapset;
- Site *site;
+ FCELL *cell1, *cell4b,*cell5;
+ FCELL *cell9, *cell10,*cell11;
+ DCELL *cell2, *cell3, *cell4, *cell4a,*cell12;
+ int fd1, fd2, fd3, fd4, fd4a, fd4b, fd5, row, row_rev;
+ int fd9, fd10, fd11, fd12;
+ int l,j;
+ int nn,cc,ii,dd;
+ double unitconv=0.0000002; /* mm/hr to m/s */
+ char *mapset;
+ Site *site;
- npoints = 0;
- npoints_alloc = 0;
+ npoints = 0;
+ npoints_alloc = 0;
-/* put some warning messages about diff. in resol., etc. */
+ /* put some warning messages about diff. in resol., etc. */
+ if (sfile != NULL) {
+ /* open ascii file for ascii output of gamma */
+ fw=fopen("simwe_data.txt","w");
- if (sfile != NULL) {
- fw=fopen("simwe_data.txt","w"); /* open ascii file for ascii output of gamma */
+ mapset = G_find_sites (sfile, "");
+ if (mapset == NULL)
+ G_fatal_error (_("File [%s] not found"), sfile);
- mapset = G_find_sites (sfile, "");
- if (mapset == NULL)
- G_fatal_error (_("File [%s] not found"), sfile);
+ if ((fdsfile = G_fopen_sites_old (sfile, mapset)) == NULL)
+ G_fatal_error (_("Unable to open file [%s]"), sfile);
- if ((fdsfile = G_fopen_sites_old (sfile, mapset)) == NULL)
- G_fatal_error (_("Unable to open file [%s]"), sfile);
+ if (G_site_describe (fdsfile, &nn, &cc, &ii, &dd)!=0)
+ G_fatal_error (_("Failed to guess file format"));
- if (G_site_describe (fdsfile, &nn, &cc, &ii, &dd)!=0)
- G_fatal_error (_("Failed to guess file format"));
+ site = G_site_new_struct (cc, nn, ii, dd);
+ G_message (_("Reading sites map (%s) ..."), sfile);
- site = G_site_new_struct (cc, nn, ii, dd);
- G_message (_("Reading sites map (%s) ..."), sfile);
+/*
+ if (dd==0)
+ {
+ fprintf(stderr,"\n");
+ G_warning("I'm finding records that do not have
+ a floating point attributes (fields prefixed with '%').");
+ }
+*/
-/* if (dd==0)
- {
- fprintf(stderr,"\n");
- G_warning("I'm finding records that do not have a floating point attributes (fields prefixed with '%').");
- } */
-
- while (G_site_get(fdsfile, site) >= 0)
- {
- if (npoints_alloc <= npoints)
- {
- npoints_alloc += 128;
- points = (struct Point*) G_realloc(points, npoints_alloc * sizeof (struct Point));
- }
-
- points[npoints].east = site->east * conv;
- points[npoints].north = site->north * conv;
- points[npoints].z1 = 0.; /*site->dbl_att[0];*/
-
-/*printf("\n%f %f",points[npoints].east/conv,points[npoints].north/conv);*/
-if((points[npoints].east/conv <= cellhd.east && points[npoints].east/conv >= cellhd.west) &&
- (points[npoints].north/conv <= cellhd.north && points[npoints].north/conv >= cellhd.south))
- npoints++;
- }
-
-G_sites_close(fdsfile);
+ while (G_site_get(fdsfile, site) >= 0)
+ {
+ if (npoints_alloc <= npoints)
+ {
+ npoints_alloc += 128;
+ points = (struct Point*) G_realloc(points, npoints_alloc * sizeof (struct Point));
+ }
+ points[npoints].east = site->east * conv;
+ points[npoints].north = site->north * conv;
+ points[npoints].z1 = 0.; /*site->dbl_att[0];*/
+ /*printf("\n%f %f",points[npoints].east/conv,points[npoints].north/conv);*/
+ if((points[npoints].east/conv <= cellhd.east && points[npoints].east/conv >= cellhd.west) && (points[npoints].north/conv <= cellhd.north && points[npoints].north/conv >= cellhd.south))
+ npoints++;
+ }
+ G_sites_close(fdsfile);
}
+ /* Allocate raster buffers */
cell1=G_allocate_f_raster_buf();
cell2=G_allocate_d_raster_buf();
cell3=G_allocate_d_raster_buf();
- if(rain != NULL)
- cell4=G_allocate_d_raster_buf();
+ if(rain != NULL)
+ cell4=G_allocate_d_raster_buf();
- if(infil != NULL)
- cell4a=G_allocate_d_raster_buf();
+ if(infil != NULL)
+ cell4a=G_allocate_d_raster_buf();
- if(traps != NULL)
- cell4b=G_allocate_f_raster_buf();
- cell5=G_allocate_f_raster_buf();
+ if(traps != NULL)
+ cell4b=G_allocate_f_raster_buf();
+ cell5=G_allocate_f_raster_buf();
- if(detin!=NULL)
- cell9=G_allocate_f_raster_buf();
+ if(detin!=NULL)
+ cell9=G_allocate_f_raster_buf();
- if(tranin!=NULL)
- cell10=G_allocate_f_raster_buf();
+ if(tranin!=NULL)
+ cell10=G_allocate_f_raster_buf();
- if(tauin!=NULL)
- cell11=G_allocate_f_raster_buf();
+ if(tauin!=NULL)
+ cell11=G_allocate_f_raster_buf();
- if(wdepth!=NULL)
- cell12=G_allocate_d_raster_buf();
+ if(wdepth!=NULL)
+ cell12=G_allocate_d_raster_buf();
- zz = (float **)G_malloc (sizeof(float *)*(my));
- v1 = (double **)G_malloc (sizeof(double *)*(my));
- v2 = (double **)G_malloc (sizeof(double *)*(my));
+ /* Allocate some double dimension arrays for each input*/
+ /* with length of matrix Y */
+ zz = (float **)G_malloc (sizeof(float *)*(my));
+ v1 = (double **)G_malloc (sizeof(double *)*(my));
+ v2 = (double **)G_malloc (sizeof(double *)*(my));
- if(rain != NULL)
- si = (double **)G_malloc (sizeof(double *)*(my));
+ if(rain != NULL||rain_val >= 0.0)
+ si = (double **)G_malloc (sizeof(double *)*(my));
- if(infil != NULL)
- inf = (double **)G_malloc (sizeof(double *)*(my));
+ if(infil != NULL||infil_val >= 0.0)
+ inf = (double **)G_malloc (sizeof(double *)*(my));
- if(traps != NULL)
- trap = (float **)G_malloc (sizeof(float *)*(my));
- cchez = (float **)G_malloc (sizeof(float *)*(my));
+ if(traps != NULL)
+ trap = (float **)G_malloc (sizeof(float *)*(my));
+ cchez = (float **)G_malloc (sizeof(float *)*(my));
- if(detin!=NULL)
- dc = (float **)G_malloc (sizeof(float *)*(my));
+ if(detin!=NULL)
+ dc = (float **)G_malloc (sizeof(float *)*(my));
- if(tranin!=NULL)
- ct = (float **)G_malloc (sizeof(float *)*(my));
+ if(tranin!=NULL)
+ ct = (float **)G_malloc (sizeof(float *)*(my));
- if(tauin!=NULL)
- tau = (float **)G_malloc (sizeof(float *)*(my));
+ if(tauin!=NULL)
+ tau = (float **)G_malloc (sizeof(float *)*(my));
+
+ if(wdepth!=NULL)
+ gama = (double **)G_malloc (sizeof(double *)*(my));
- if(wdepth!=NULL)
- gama = (double **)G_malloc (sizeof(double *)*(my));
+ for(l=0;l<my;l++)
+ {
+ /*for each my, allocate a second dimension in array
+ * for each input with length of matrix X*/
+ zz[l] = (float*)G_malloc (sizeof(float)*(mx));
+ v1[l] = (double*)G_malloc (sizeof(double)*(mx));
+ v2[l] = (double*)G_malloc (sizeof(double)*(mx));
- for(l=0;l<my;l++)
- {
- zz[l] = (float*)G_malloc (sizeof(float)*(mx));
- v1[l] = (double*)G_malloc (sizeof(double)*(mx));
- v2[l] = (double*)G_malloc (sizeof(double)*(mx));
+ if(rain != NULL||rain_val >= 0.0)
+ si[l] = (double*)G_malloc (sizeof(double)*(mx));
- if(rain != NULL)
- si[l] = (double*)G_malloc (sizeof(double)*(mx));
+ if(infil != NULL||infil_val >= 0.0)
+ inf[l] = (double*)G_malloc (sizeof(double)*(mx));
- if(infil != NULL)
- inf[l] = (double*)G_malloc (sizeof(double)*(mx));
+ if(traps != NULL)
+ trap[l] = (float*)G_malloc (sizeof(float)*(mx));
+ cchez[l] = (float*)G_malloc (sizeof(float)*(mx));
- if(traps != NULL)
- trap[l] = (float*)G_malloc (sizeof(float)*(mx));
- cchez[l] = (float*)G_malloc (sizeof(float)*(mx));
+ if(detin!=NULL)
+ dc[l] = (float*)G_malloc (sizeof(float)*(mx));
- if(detin!=NULL)
- dc[l] = (float*)G_malloc (sizeof(float)*(mx));
+ if(tranin!=NULL)
+ ct[l] = (float*)G_malloc (sizeof(float)*(mx));
- if(tranin!=NULL)
- ct[l] = (float*)G_malloc (sizeof(float)*(mx));
+ if(tauin!=NULL)
+ tau[l] = (float*)G_malloc (sizeof(float)*(mx));
- if(tauin!=NULL)
- tau[l] = (float*)G_malloc (sizeof(float)*(mx));
+ if(wdepth!=NULL)
+ gama[l] = (double*)G_malloc (sizeof(double)*(mx));
+ }
- if(wdepth!=NULL)
- gama[l] = (double*)G_malloc (sizeof(double)*(mx));
- }
+ G_message (_("Running MAY 10 version, started modifications on 20080211"));
- G_message (_("Running MAY 10 version"));
+ /* Check if data available in mapsets
+ * if found, then open the files */
+ if((mapset=G_find_cell(elevin,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), elevin);
- if((mapset=G_find_cell(elevin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), elevin);
+ fd1 = G_open_cell_old(elevin,mapset);
- fd1 = G_open_cell_old(elevin,mapset);
+ /* TO REPLACE BY INTERNAL PROCESSING of dx, dy from Elevin*/
+ if((mapset=G_find_cell(dxin,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), dxin);
- if((mapset=G_find_cell(dxin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), dxin);
+ fd2 = G_open_cell_old(dxin,mapset);
- fd2 = G_open_cell_old(dxin,mapset);
+ if((mapset=G_find_cell(dyin,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), dyin);
- if((mapset=G_find_cell(dyin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), dyin);
+ fd3 = G_open_cell_old(dyin,mapset);
+ /* END OF REPLACEMENT */
- fd3 = G_open_cell_old(dyin,mapset);
+ /* Rendered Mannings n input map optional to run!*/
+ /* Careful! (Yann, 20080212)*/
+ if(manin != NULL){
+ if((mapset=G_find_cell(manin,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), manin);
+ fd5 = G_open_cell_old(manin,mapset);
+ }
- if((mapset=G_find_cell(manin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), manin);
+ /* Rendered Rainfall input map optional to run!*/
+ /* Careful! (Yann, 20080212)*/
+ if(rain != NULL){
+ if((mapset=G_find_cell(rain,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), rain);
+ fd4 = G_open_cell_old(rain,mapset);
+ }
- fd5 = G_open_cell_old(manin,mapset);
+ if(infil != NULL){
+ if((mapset=G_find_cell(infil,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), infil);
+ fd4a = G_open_cell_old(infil,mapset);
+ }
- if(rain != NULL){
- if((mapset=G_find_cell(rain,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), rain);
- fd4 = G_open_cell_old(rain,mapset);
- }
+ if(traps != NULL){
+ if((mapset=G_find_cell(traps,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), traps);
- if(infil != NULL){
- if((mapset=G_find_cell(infil,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), infil);
+ fd4b = G_open_cell_old(traps,mapset);
+ }
- fd4a = G_open_cell_old(infil,mapset);
- }
+ if(detin != NULL){
+ if((mapset=G_find_cell(detin,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), detin);
- if(traps != NULL){
- if((mapset=G_find_cell(traps,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), traps);
+ fd9 = G_open_cell_old(detin,mapset);
+ }
- fd4b = G_open_cell_old(traps,mapset);
- }
+ if(tranin != NULL){
+ if((mapset=G_find_cell(tranin,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), tranin);
- if(detin != NULL){
- if((mapset=G_find_cell(detin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), detin);
+ fd10 = G_open_cell_old(tranin,mapset);
+ }
- fd9 = G_open_cell_old(detin,mapset);
- }
+ if(tauin != NULL){
+ if((mapset=G_find_cell(tauin,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), tauin);
- if(tranin != NULL){
- if((mapset=G_find_cell(tranin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), tranin);
+ fd11 = G_open_cell_old(tauin,mapset);
+ }
- fd10 = G_open_cell_old(tranin,mapset);
- }
+ if(wdepth != NULL){
+ if((mapset=G_find_cell(wdepth,""))==NULL)
+ G_fatal_error(_("Raster map <%s> not found"), wdepth);
- if(tauin != NULL){
- if((mapset=G_find_cell(tauin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), tauin);
+ fd12 = G_open_cell_old(wdepth,mapset);
+ }
- fd11 = G_open_cell_old(tauin,mapset);
- }
+ for (row=0; row<my; row++)
+ {
+ G_get_f_raster_row(fd1,cell1,row);
+ G_get_d_raster_row(fd2,cell2,row);
+ G_get_d_raster_row(fd3,cell3,row);
- if(wdepth != NULL){
- if((mapset=G_find_cell(wdepth,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), wdepth);
+ if(manin != NULL)
+ G_get_f_raster_row(fd5,cell5,row);
- fd12 = G_open_cell_old(wdepth,mapset);
- }
+ if(rain != NULL)
+ G_get_d_raster_row(fd4,cell4,row);
- for (row=0; row<my; row++)
- {
- G_get_f_raster_row(fd1,cell1,row);
- G_get_d_raster_row(fd2,cell2,row);
- G_get_d_raster_row(fd3,cell3,row);
+ if(infil != NULL)
+ G_get_d_raster_row(fd4a,cell4a,row);
- if(rain != NULL)
- G_get_d_raster_row(fd4,cell4,row);
+ if(traps != NULL)
+ G_get_f_raster_row(fd4b,cell4b,row);
- if(infil != NULL)
- G_get_d_raster_row(fd4a,cell4a,row);
+ if(detin != NULL)
+ G_get_f_raster_row(fd9,cell9,row);
- if(traps != NULL)
- G_get_f_raster_row(fd4b,cell4b,row);
- G_get_f_raster_row(fd5,cell5,row);
+ if(tranin != NULL)
+ G_get_f_raster_row(fd10,cell10,row);
- if(detin != NULL)
- G_get_f_raster_row(fd9,cell9,row);
+ if(tauin != NULL)
+ G_get_f_raster_row(fd11,cell11,row);
- if(tranin != NULL)
- G_get_f_raster_row(fd10,cell10,row);
+ if(wdepth != NULL)
+ G_get_d_raster_row(fd12,cell12,row);
- if(tauin != NULL)
- G_get_f_raster_row(fd11,cell11,row);
+ for (j=0; j<mx; j++)
+ {
+ row_rev = my - row - 1;
- if(wdepth != NULL)
- G_get_d_raster_row(fd12,cell12,row);
+ if(!G_is_f_null_value(cell1+j))
+ zz[row_rev][j] = (float ) (conv * cell1[j]);
+ else
+ zz[row_rev][j] = UNDEF;
- for (j=0; j<mx; j++)
- {
- row_rev = my - row - 1;
+ if(!G_is_d_null_value(cell2+j))
+ v1[row_rev][j] = (double ) cell2[j];
+ else
+ v1[row_rev][j] = UNDEF;
- if(!G_is_f_null_value(cell1+j))
- zz[row_rev][j] = (float ) (conv * cell1[j]);
- else
- zz[row_rev][j] = UNDEF;
+ if(!G_is_d_null_value(cell3+j))
+ v2[row_rev][j] = (double ) cell3[j];
+ else
+ v2[row_rev][j] = UNDEF;
- if(!G_is_d_null_value(cell2+j))
- v1[row_rev][j] = (double ) cell2[j];
- else
- v1[row_rev][j] = UNDEF;
-
- if(!G_is_d_null_value(cell3+j))
- v2[row_rev][j] = (double ) cell3[j];
- else
- v2[row_rev][j] = UNDEF;
-
- if(v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF)
- zz[row_rev][j] = UNDEF; /* undef all area if something's missing */
+ if(v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF)
+ zz[row_rev][j] = UNDEF; /* undef all area if something's missing */
-/* should be ? if(v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF || zz[row_rev][j] = UNDEF)
- * {v1[row_rev][j] == UNDEF;
- * v2[row_rev][j] == UNDEF;
- * zz[row_rev][j] = UNDEF;}
- * printout warning?
- */
- if (rain != NULL)
- {
- if(!G_is_d_null_value(cell4+j))
- si[row_rev][j] = ((double ) cell4[j]) * unitconv; /*conv mm/hr to m/s*/
- /* printf("\n INPUTrain, convert %f %f",si[row_rev][j],unitconv); */
+ /* should be ?
+ * if(v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF ||
+ * zz[row_rev][j] = UNDEF)
+ * {v1[row_rev][j] == UNDEF;
+ * v2[row_rev][j] == UNDEF;
+ * zz[row_rev][j] = UNDEF;}
+ * printout warning?
+ */
+ if (rain != NULL)
+ {
+ if(!G_is_d_null_value(cell4+j))
+ si[row_rev][j] = ((double ) cell4[j]) * unitconv; /*conv mm/hr to m/s*/
+ /*printf("\n INPUTrain, convert %f %f",si[row_rev][j],unitconv); */
- else {
- si[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
+ else {
+ si[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
- if (infil != NULL)
- {
- if(!G_is_d_null_value(cell4a+j))
- inf[row_rev][j] = (double ) cell4a[j] * unitconv; /*conv mm/hr to m/s*/
- /* printf("\n INPUT infilt, convert %f %f",inf[row_rev][j],unitconv); */
- else {
- inf[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
+ if (infil != NULL)
+ {
+ if(!G_is_d_null_value(cell4a+j))
+ inf[row_rev][j] = (double ) cell4a[j] * unitconv; /*conv mm/hr to m/s*/
+ /*printf("\nINPUT infilt,convert %f %f",inf[row_rev][j],unitconv);*/
+ else {
+ inf[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ } else { /* Added by Yann 20080216*/
+ /* If infil==NULL, then use infilval */
+ if( infil_val >= 0.0 ){
+ inf[row_rev][j]= infil_val * unitconv; /*conv mm/hr to m/s*/
+ /* printf("infil_val = %f \n",inf[row_rev][j]);*/
+ } else {
+ inf[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
- if (traps != NULL)
- {
- if(!G_is_f_null_value(cell4b+j))
- trap[row_rev][j] = (float) cell4b[j]; /* no conv, unitless */
- else {
- trap[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- }
+ if (traps != NULL)
+ {
+ if(!G_is_f_null_value(cell4b+j))
+ trap[row_rev][j] = (float) cell4b[j]; /* no conv, unitless */
+ else {
+ trap[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
+ } else { /* Added by Yann 20080213*/
+ /* If rain==NULL, then use rainval */
+ if(rain_val>0.0){
+ si[row_rev][j]= rain_val*unitconv; /* conv mm/hr to m/s */
+ /*printf("\n INPUTrainval, convert %f %f",si[row_rev][j],unitconv); */
+ } else {
+ si[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+
+ if (infil != NULL)
+ {
+ if(!G_is_d_null_value(cell4a+j))
+ inf[row_rev][j] = (double ) cell4a[j] * unitconv; /*conv mm/hr to m/s*/
+ /*printf("\nINPUT infilt,convert %f %f",inf[row_rev][j],unitconv);*/
+ else {
+ inf[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ } else { /* Added by Yann 20080216*/
+ /* If infil==NULL, then use infilval */
+ if( infil_val >= 0.0 ){
+ inf[row_rev][j]= infil_val * unitconv; /*conv mm/hr to m/s*/
+ /*printf("infil_val = %f \n",inf[row_rev][j]);*/
+ } else {
+ inf[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
+
+ if (traps != NULL)
+ {
+ if(!G_is_f_null_value(cell4b+j))
+ trap[row_rev][j] = (float) cell4b[j]; /* no conv, unitless */
+ else {
+ trap[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
+ } /* End of added by Yann 20080213*/
+ if (manin != NULL){
+ if(!G_is_f_null_value(cell5+j)){
+ cchez[row_rev][j] = (float ) cell5[j]; /* units in manual */
+ } else {
+ cchez[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ } else if(manin_val>0.0) { /* Added by Yann 20080213 */
+ cchez[row_rev][j] = (float) manin_val;
+ } else {
+ G_fatal_error(_("Raster map <%s> not found, and manin_val undefined, choose one to be allowed to process"), manin);
+ }
+ if (detin != NULL)
+ {
+ if(!G_is_f_null_value(cell9+j))
+ dc[row_rev][j] = (float ) cell9[j]; /*units in manual*/
+ else {
+ dc[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
- if(!G_is_f_null_value(cell5+j))
- cchez[row_rev][j] = (float ) cell5[j]; /* units in manual */
- else {
- cchez[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
+ if (tranin != NULL)
+ {
+ if(!G_is_f_null_value(cell10+j))
+ ct[row_rev][j] = (float ) cell10[j]; /*units in manual*/
+ else {
+ ct[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
- if (detin != NULL)
- {
- if(!G_is_f_null_value(cell9+j))
- dc[row_rev][j] = (float ) cell9[j]; /*units in manual*/
- else {
- dc[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
+ if (tauin != NULL)
+ {
+ if(!G_is_f_null_value(cell11+j))
+ tau[row_rev][j] = (float ) cell11[j]; /*units in manual*/
+ else {
+ tau[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
- if (tranin != NULL)
- {
- if(!G_is_f_null_value(cell10+j))
- ct[row_rev][j] = (float ) cell10[j]; /*units in manual*/
- else {
- ct[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
+ if (wdepth != NULL)
+ {
+ if(!G_is_d_null_value(cell12+j))
+ gama[row_rev][j] = (double ) cell12[j]; /*units in manual*/
+ else {
+ gama[row_rev][j] = UNDEF;
+ zz[row_rev][j] = UNDEF;
+ }
+ }
}
- }
+ }
+ G_close_cell(fd1);
+ G_close_cell(fd2);
+ G_close_cell(fd3);
- if (tauin != NULL)
- {
- if(!G_is_f_null_value(cell11+j))
- tau[row_rev][j] = (float ) cell11[j]; /*units in manual*/
- else {
- tau[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
+ if(rain != NULL)
+ G_close_cell(fd4);
- if (wdepth != NULL)
- {
- if(!G_is_d_null_value(cell12+j))
- gama[row_rev][j] = (double ) cell12[j]; /*units in manual*/
- else {
- gama[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- }
- }
- G_close_cell(fd1);
- G_close_cell(fd2);
- G_close_cell(fd3);
+ if(infil != NULL)
+ G_close_cell(fd4a);
- if(rain != NULL)
- G_close_cell(fd4);
+ if(traps != NULL)
+ G_close_cell(fd4b);
+ /* Maybe a conditional to manin!=NULL here ! */
+ G_close_cell(fd5);
+ /****************/
- if(infil != NULL)
- G_close_cell(fd4a);
+ if(detin != NULL)
+ G_close_cell(fd9);
- if(traps != NULL)
- G_close_cell(fd4b);
- G_close_cell(fd5);
+ if(tranin != NULL)
+ G_close_cell(fd10);
- if(detin != NULL)
- G_close_cell(fd9);
+ if(tauin != NULL)
+ G_close_cell(fd11);
- if(tranin != NULL)
- G_close_cell(fd10);
+ if(wdepth != NULL)
+ G_close_cell(fd12);
- if(tauin != NULL)
- G_close_cell(fd11);
-
- if(wdepth != NULL)
- G_close_cell(fd12);
-
- return 1;
+ return 1;
}
/* data preparations, sigma, shear, etc. */
int grad_check (void)
{
- int k,l,i,j;
- double zx,zy,zd2,zd4,sinsl;
-
- double cc,cmul2;
- double sheer;
- double vsum = 0.;
- double vmax = 0.;
- double chsum = 0.;
- double zmin = 1.e12;
- double zmax = -1.e12;
- double zd2min = 1.e12;
- double zd2max = -1.e12;
- double smin = 1.e12;
- double smax = -1.e12;
- double infmin = 1.e12;
- double infmax = -1.e12;
- double sigmax = -1.e12;
- double cchezmax = -1.e12;
- double rhow = 1000.;
- double gacc = 9.81;
- double hh = 1.;
- double deltaw = 1.e12;
+ int k,l,i,j;
+ double zx,zy,zd2,zd4,sinsl;
+ double cc,cmul2;
+ double sheer;
+ double vsum = 0.;
+ double vmax = 0.;
+ double chsum = 0.;
+ double zmin = 1.e12;
+ double zmax = -1.e12;
+ double zd2min = 1.e12;
+ double zd2max = -1.e12;
+ double smin = 1.e12;
+ double smax = -1.e12;
+ double infmin = 1.e12;
+ double infmax = -1.e12;
+ double sigmax = -1.e12;
+ double cchezmax = -1.e12;
+ double rhow = 1000.;
+ double gacc = 9.81;
+ double hh = 1.;
+ double deltaw = 1.e12;
sisum = 0.;
infsum = 0.;
cmul2 = rhow * gacc;
-/* mandatory alloc. - should be moved to main.c*/
+ /* mandatory alloc. - should be moved to main.c*/
slope = (double **)G_malloc (sizeof(double *)*(my));
- for(l=0;l<my;l++)
- slope[l] = (double*)G_malloc (sizeof(double)*(mx));
+ for(l=0;l<my;l++)
+ slope[l] = (double*)G_malloc (sizeof(double)*(mx));
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- slope[j][i] = 0.;
- }
-/*** */
+ for (j = 0; j < my; j++)
+ {
+ for (i = 0; i < mx; i++)
+ slope[j][i] = 0.;
+ }
+ /*** */
- for (k = 0; k < my; k++) {
- for (l = 0; l < mx; l++) {
-
- if (zz[k][l] != UNDEF) {
- zx = v1[k][l];
- zy = v2[k][l];
- zd2 = zx * zx + zy * zy;
- sinsl = sqrt(zd2) / sqrt(zd2 + 1); /* sin(terrain slope) */
-/* Computing MIN */
- zd2 = sqrt(zd2);
- zd2min = amin1(zd2min,zd2);
-/* Computing MAX */
- zd2max = amax1(zd2max,zd2);
- zd4 = sqrt(zd2); /* ^.25 */
-
- if (cchez[k][l] != 0.) {
-
- cchez[k][l] = 1. / cchez[k][l]; /* 1/n */
-
- } else {
- G_fatal_error (_("Zero value in Mannings n"));
- }
-
- if (zd2 == 0.) {
- v1[k][l] = 0.;
- v2[k][l] = 0.;
- slope[k][l] = 0.;
- } else {
- if (wdepth != NULL)
- hh = pow(gama[k][l],2./3.);
- v1[k][l] = (double)hh * cchez[k][l] * zx / zd4; /* hh = 1 if there is no water depth input */
- v2[k][l] = (double)hh * cchez[k][l] * zy / zd4;
- slope[k][l] = sqrt(v1[k][l] * v1[k][l] + v2[k][l] * v2[k][l]);
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ if (zz[k][l] != UNDEF) {
+ zx = v1[k][l];
+ zy = v2[k][l];
+ zd2 = zx * zx + zy * zy;
+ sinsl = sqrt(zd2) / sqrt(zd2 + 1); /* sin(terrain slope) */
+ /* Computing MIN */
+ zd2 = sqrt(zd2);
+ zd2min = amin1(zd2min,zd2);
+ /* Computing MAX */
+ zd2max = amax1(zd2max,zd2);
+ zd4 = sqrt(zd2); /* ^.25 */
+ if (cchez[k][l] != 0.) {
+ cchez[k][l] = 1. / cchez[k][l]; /* 1/n */
+ } else {
+ G_fatal_error (_("Zero value in Mannings n"));
+ }
+ if (zd2 == 0.) {
+ v1[k][l] = 0.;
+ v2[k][l] = 0.;
+ slope[k][l] = 0.;
+ } else {
+ if (wdepth != NULL)
+ hh = pow(gama[k][l],2./3.);
+ /* hh = 1 if there is no water depth input */
+ v1[k][l] = (double)hh * cchez[k][l] * zx / zd4;
+ v2[k][l] = (double)hh * cchez[k][l] * zy / zd4;
+ slope[k][l] = sqrt(v1[k][l] * v1[k][l] + v2[k][l] * v2[k][l]);
+ }
+ if (wdepth != NULL) {
+ sheer = (double) (cmul2 * gama[k][l] * sinsl); /* shear stress */
+ /* if critical shear stress >= shear then all zero */
+ if ((sheer <= tau[k][l]) || (ct[k][l] == 0.)) {
+ si[k][l] = 0.;
+ sigma[k][l] = 0.;
+ } else {
+ si[k][l] = (double) (dc[k][l] * (sheer - tau[k][l]));
+ sigma[k][l] = (double) (dc[k][l]/ct[k][l])*(sheer-tau[k][l])/(pow(sheer,1.5)); /* rill erosion=1.5, sheet = 1.1 */
+ }
+ }
+ sisum += si[k][l];
+ smin = amin1(smin,si[k][l]);
+ smax = amax1(smax,si[k][l]);
+ if (inf != NULL) {
+ infsum += inf[k][l];
+ infmin = amin1(infmin,inf[k][l]);
+ infmax = amax1(infmax,inf[k][l]);
+ }
+ vmax = amax1(vmax,slope[k][l]);
+ vsum += slope[k][l];
+ chsum += cchez[k][l];
+ zmin = amin1(zmin,(double)zz[k][l]);
+ zmax = amax1(zmax,(double)zz[k][l]); /* not clear were needed */
+ if (wdepth != NULL)
+ sigmax = amax1(sigmax,sigma[k][l]);
+ cchezmax = amax1(cchezmax,cchez[k][l]);
+ /* saved sqrt(sinsl)*cchez to cchez array for output*/
+ cchez[k][l] *= sqrt(sinsl);
+ } /* DEFined area */
}
-
- if (wdepth != NULL) {
- sheer = (double) (cmul2 * gama[k][l] * sinsl); /* shear stress */
-
- if ((sheer <= tau[k][l]) || (ct[k][l] == 0.)) { /* if critical shear stress >= shear then all zero */
- si[k][l] = 0.;
- sigma[k][l] = 0.;
- } else
- {
- si[k][l] = (double) (dc[k][l] * (sheer - tau[k][l]));
- sigma[k][l] = (double) (dc[k][l]/ct[k][l])*(sheer-tau[k][l])/(pow(sheer,1.5)); /* rill erosion=1.5, sheet = 1.1 */
- }
}
+ if (inf != NULL && smax < infmax)
+ G_warning (_("Infiltration exceeds the rainfall rate everywhere! No overland flow."));
- sisum += si[k][l];
- smin = amin1(smin,si[k][l]);
- smax = amax1(smax,si[k][l]);
+ cc = (double) mx * my;
+ si0 = sisum / cc;
+ vmean = vsum / cc;
+ chmean = chsum / cc;
- if (inf != NULL) {
- infsum += inf[k][l];
- infmin = amin1(infmin,inf[k][l]);
- infmax = amax1(infmax,inf[k][l]);
- }
- vmax = amax1(vmax,slope[k][l]);
- vsum += slope[k][l];
- chsum += cchez[k][l];
+ if (inf != NULL) infmean = infsum / cc;
- zmin = amin1(zmin,(double)zz[k][l]);
- zmax = amax1(zmax,(double)zz[k][l]); /* not clear were needed */
+ if (wdepth != NULL) deltaw = 0.8/(sigmax * vmax); /*time step for sediment*/
+ deltap = 0.25 * sqrt(stepx * stepy)/vmean; /*time step for water*/
- if (wdepth != NULL)
- sigmax = amax1(sigmax,sigma[k][l]);
- cchezmax = amax1(cchezmax,cchez[k][l]);
+ if(deltaw > deltap)
+ timec = 4.;
+ else
+ timec = 1.25;
- cchez[k][l] *= sqrt(sinsl); /* saved sqrt(sinsl)*cchez to cchez array for output*/
- } /* DEFined area */
- }
- }
+ miter = (int)(timesec / (deltap * timec)); /* number of iterations = number of cells to pass*/
+ iterout = (int)(iterout / (deltap * timec)); /* number of cells to pass for time series output */
+
+ fprintf (stderr, "\n");
+ G_message (_("Min elevation \t= %.2f m\nMax elevation \t= %.2f m\n"),zmin,zmax);
+ G_message (_("Mean Rainfall \t= %f m/s\nMean flow velocity \t= %f m/s\n"),si0,vmean);
+ G_message (_("Mean Mannings \t= %f\n"),1.0/chmean);
- if (inf != NULL && smax < infmax)
- G_warning (_("Infiltration exceeds the rainfall rate everywhere! No overland flow."));
+ deltap = amin1(deltap,deltaw);
- cc = (double) mx * my;
- si0 = sisum / cc;
- vmean = vsum / cc;
- chmean = chsum / cc;
+ G_message (_("Number of iterations \t= %d cells\n"),miter);
+ G_message (_("Time step \t= %.2f s\n"),deltap);
+ if(wdepth != NULL){
+ G_message (_("Sigmax \t= %f\nMax velocity \t= %f m/s\n"),sigmax,vmax);
+ G_message (_("Time step used \t= %.2f s\n"), deltaw);
+ }
+ /* if (wdepth != NULL) deltap = 0.1;
+ * deltap for sediment is ar. average deltap and deltaw */
+ /* if (wdepth != NULL) deltap = (deltaw+deltap)/2.;
+ * deltap for sediment is ar. average deltap and deltaw */
- if (inf != NULL) infmean = infsum / cc;
- if (wdepth != NULL) deltaw = 0.8/(sigmax * vmax);
- deltap = 0.25 * sqrt(stepx * stepy)/vmean;
-
- if(deltaw > deltap)
- timec = 4.;
- else
- timec = 1.25;
-
- fprintf (stderr, "\n");
- G_message ("zmin,zmax %f %f", zmin, zmax);
- G_message ("simean,vmean,chmean,deltap,deltaw %f %f %f %f %f",
- si0, vmean, chmean, deltap, deltaw);
- G_message ("MITER, timec %d %f", miter, timec);
-
- deltap = amin1(deltap,deltaw);
-
- if (wdepth != NULL)
- G_message ("sigmax,vmax,deltapused %f %f %f", sigmax, vmax, deltaw);
-
-/* if (wdepth != NULL) deltap = 0.1; deltap for sediment is ar. average deltap and deltaw */
-/* if (wdepth != NULL) deltap = (deltaw+deltap)/2.; deltap for sediment is ar. average deltap and deltaw */
-
- miter = (int)(timesec / (deltap * timec)); /* number of iterations */
- iterout = (int)(iterout / (deltap * timec)); /* iterations for time series */
-
-/*! For each cell (k,l) compute the length s=(v1,v2) of the path
- * that the particle will travel per one time step
- * \f$ s(k,l)=v(k,l)*dt \f$, [m]=[m/s]*[s]
- * give warning if there is a cell that will lead to path longer than 2 cells
- *
- * if running erosion, compute sediment transport capacity for each cell si(k,l)
- * \f$
- * T({\bf r})=K_t({\bf r}) \bigl[\tau({\bf r})\bigr]^p
- * =K_t({\bf r}) \bigl[\rho_w\, g h({\bf r}) \sin \beta ({\bf r}) \bigr]^p
- * \f$
- * [kg/ms]=...
-*/
- for (k = 0; k < my; k++) {
- for (l = 0; l < mx; l++) {
-
- if (zz[k][l] != UNDEF) {
- v1[k][l] *= deltap;
- v2[k][l] *= deltap;
-
-/* if(v1[k][l]*v1[k][l]+v2[k][l]*v2[k][l] > cellsize, warning, napocitaj
-* ak viac ako 10%a*/
-
- if (inf!=NULL) inf[k][l] *= timesec; /* THIS IS CORRECT SOLUTION currently commented out*/
-
- if(wdepth != NULL) gama[k][l] = 0.;
-
- if(et!=NULL) {
- if(sigma[k][l] == 0. || slope[k][l] == 0.)
- si[k][l] = 0.;
- else
- si[k][l] = si[k][l] / (slope[k][l] * sigma[k][l]); /* temp for transp. cap. erod */
+ /*! For each cell (k,l) compute the length s=(v1,v2) of the path
+ * that the particle will travel per one time step
+ * \f$ s(k,l)=v(k,l)*dt \f$, [m]=[m/s]*[s]
+ * give warning if there is a cell that will lead to path longer than 2 cells
+ *
+ * if running erosion, compute sediment transport capacity for each cell si(k,l)
+ * \f$
+ * T({\bf r})=K_t({\bf r}) \bigl[\tau({\bf r})\bigr]^p
+ * =K_t({\bf r}) \bigl[\rho_w\, g h({\bf r}) \sin \beta ({\bf r}) \bigr]^p
+ * \f$
+ * [kg/ms]=...
+ */
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ if (zz[k][l] != UNDEF) {
+ v1[k][l] *= deltap;
+ v2[k][l] *= deltap;
+ /*if(v1[k][l]*v1[k][l]+v2[k][l]*v2[k][l] > cellsize, warning, napocitaj
+ *ak viac ako 10%a*/
+ /* THIS IS CORRECT SOLUTION currently commented out*/
+ if (inf!=NULL) inf[k][l] *= timesec;
+ if(wdepth != NULL) gama[k][l] = 0.;
+ if(et!=NULL) {
+ if(sigma[k][l] == 0. || slope[k][l] == 0.)
+ si[k][l] = 0.;
+ else
+ /* temp for transp. cap. erod */
+ si[k][l] = si[k][l] / (slope[k][l] * sigma[k][l]);
+ }
+ } /* DEFined area */
}
- } /* DEFined area */
- }
- }
+ }
-/*! compute transport capacity limted erosion/deposition et
-* as a divergence of sediment transport capacity
-* \f$
- D_T({\bf r})= \nabla\cdot {\bf T}({\bf r})
-* \f$
-*/
+ /*! compute transport capacity limted erosion/deposition et
+ * as a divergence of sediment transport capacity
+ * \f$
+ D_T({\bf r})= \nabla\cdot {\bf T}({\bf r})
+ * \f$
+ */
if(et!=NULL) {
- erod(si); /* compute divergence of t.capc */
-
- if (output_et() != 1)
- G_fatal_error (_("Unable to write et file"));
+ erod(si); /* compute divergence of t.capc */
+ if (output_et() != 1)
+ G_fatal_error (_("Unable to write et file"));
}
-/*! compute the inversion operator and store it in sigma - note that after this
-* sigma does not store the first order reaction coefficient but the operator
-* WRITE the equation here
-*/
+ /*! compute the inversion operator and store it in sigma - note that after this
+ * sigma does not store the first order reaction coefficient but the operator
+ * WRITE the equation here
+ */
if(wdepth != NULL) {
- for (k = 0; k < my; k++) {
- for (l = 0; l < mx; l++) {
- if (zz[k][l] != UNDEF) {
-
- if(et!=NULL) si[k][l] = si[k][l] * slope[k][l] * sigma[k][l]; /* get back from temp */
-
- if( sigma[k][l] != 0.)
-/* rate of weight loss - w=w*sigma , vaha prechadzky po n-krokoch je sigma^n */
- sigma[k][l] = exp(-sigma[k][l] * deltap * slope[k][l]); /* not clear what's here :-\ */
-/* if(sigma[k][l]<0.5) warning, napocitaj, ak vacsie ako 50% skonci, zmensi deltap)*/
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ if (zz[k][l] != UNDEF) {
+ /* get back from temp */
+ if(et!=NULL) si[k][l] = si[k][l] * slope[k][l] * sigma[k][l];
+ if( sigma[k][l] != 0.)
+ /* rate of weight loss - w=w*sigma ,
+ * vaha prechadzky po n-krokoch je sigma^n */
+ /* not clear what's here :-\ */
+ sigma[k][l] = exp(-sigma[k][l] * deltap * slope[k][l]);
+ /* if(sigma[k][l]<0.5) warning, napocitaj,
+ * ak vacsie ako 50% skonci, zmensi deltap)*/
+ }
+ } /*DEFined area */
}
- } /*DEFined area */
- }
}
- return 1;
+ return 1;
}
double amax1 (double arg1, double arg2)
{
- double res;
+ double res;
- if (arg1>=arg2) {
- res = arg1;
- }
- else {
- res = arg2;
- }
+ if (arg1>=arg2) {
+ res = arg1;
+ } else {
+ res = arg2;
+ }
- return res;
+ return res;
}
double amin1 (double arg1, double arg2)
{
- double res;
+ double res;
- if (arg1<=arg2) {
- res = arg1;
- }
- else {
- res = arg2;
- }
+ if (arg1<=arg2) {
+ res = arg1;
+ } else {
+ res = arg2;
+ }
- return res;
+ return res;
}
int min (int arg1, int arg2)
{
- int res;
+ int res;
- if (arg1 <= arg2)
- {
- res = arg1;
- }
- else
- {
- res = arg2;
- }
+ if (arg1 <= arg2)
+ {
+ res = arg1;
+ } else {
+ res = arg2;
+ }
- return res;
+ return res;
}
int max (int arg1, int arg2)
{
- int res;
+ int res;
- if (arg1>=arg2) {
- res = arg1;
- }
- else {
- res = arg2;
- }
+ if (arg1>=arg2) {
+ res = arg1;
+ } else {
+ res = arg2;
+ }
- return res;
+ return res;
}
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/output.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/output.c 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/output.c 2008-02-18 20:52:12 UTC (rev 30240)
@@ -489,6 +489,9 @@
sprintf (hist.datsrc_1, "input files: %s %s %s", elevin,dxin,dyin);
sprintf (hist.datsrc_2, "input files: %s %s %s", rain,infil,manin);
hist.edlinecnt = 4;
+
+ G_command_history(&hist);
+
if (ts == 1) G_write_history (depth0, &hist);
else
G_write_history (depth, &hist);
@@ -508,7 +511,7 @@
/* fprintf (stdout, "\n history initiated\n");
fflush(stdout);*/
- sprintf(hist.edhist[0],"init.walk=%d, maxwalk=%d, remaining walkers=%d",nwalk,maxwa,nwalka);
+ sprintf(hist.edhist[0],"init.walkers=%d, maxwalk=%d, rem. walkers=%d",nwalk,maxwa,nwalka);
sprintf(hist.edhist[1],"duration (sec.)=%d, time-serie iteration=%d",timesec,tt);
sprintf(hist.edhist[2],"written walkers=%d, deltap=%f, mean vel.=%f",
@@ -519,6 +522,9 @@
sprintf (hist.datsrc_1, "input files: %s %s %s", elevin,dxin,dyin);
sprintf (hist.datsrc_2, "input files: %s %s %s", rain,infil,manin);
hist.edlinecnt = 4;
+
+ G_command_history(&hist);
+
if (ts == 1) G_write_history (disch0, &hist);
else
G_write_history (disch, &hist);
@@ -550,6 +556,9 @@
sprintf (hist.datsrc_2, "input files: %s %s %s %s", manin,detin,tranin,tauin);
hist.edlinecnt = 4;
+
+ G_command_history(&hist);
+
if (ts == 1) G_write_history (flux0, &hist);
else
G_write_history (flux, &hist);
@@ -704,7 +713,6 @@
/* } */
}
-
return 1;
}
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h 2008-02-18 17:34:38 UTC (rev 30239)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h 2008-02-18 20:52:12 UTC (rev 30240)
@@ -41,11 +41,15 @@
GLOBAL char *flux;
GLOBAL char *erdep;
+GLOBAL char *rainval;
+GLOBAL char *maninval;
+GLOBAL char *infilval;
+
GLOBAL struct
{
struct Option *elevin,*dxin,*dyin,*rain,*infil,*traps,*manin,*sfile,*depth,*disch,*err,
*outwalk,*nwalk,*niter,*outiter,*density,*diffc,*hmax,*halpha,*hbeta,*wdepth,
-*detin,*tranin,*tauin,*tc,*et,*conc,*flux,*erdep;
+*detin,*tranin,*tauin,*tc,*et,*conc,*flux,*erdep,*rainval,*maninval,*infilval;
} parm;
@@ -124,4 +128,10 @@
GLOBAL double timec;
GLOBAL int ts, timesec;
+GLOBAL double rain_val;
+GLOBAL double manin_val;
+GLOBAL double infil_val;
+
+GLOBAL struct History history; /* holds meta-data (title, comments,..) */
+
#endif /* __WATERGLOBS_H__ */
More information about the grass-commit
mailing list