[GRASS-SVN] r30286 - 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
Fri Feb 22 09:27:39 EST 2008
Author: neteler
Date: 2008-02-22 09:27:39 -0500 (Fri, 22 Feb 2008)
New Revision: 30286
Modified:
grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/main.c
grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/waterglobs.h
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/simlib/input.c
Log:
fixes backported
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-22 12:44:58 UTC (rev 30285)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/main.c 2008-02-22 14:27:39 UTC (rev 30286)
@@ -1,10 +1,10 @@
/****************************************************************************
*
- * MODULE: r.simwe.erosion: main program for hydrologic and sediment transport
+ * MODULE: r.sim.sediment: main program for sediment transport
* simulation (SIMWE)
*
* AUTHOR(S): L. Mitas, H. Mitasova, J. Hofierka
- * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
+ * PURPOSE: Sediment transport simulation (SIMWE)
*
* COPYRIGHT: (C) 2002 by the GRASS Development Team
*
@@ -15,7 +15,7 @@
*****************************************************************************/
/*-
- * r.simwe.erosion: main program for hydrologic and sediment transport
+ * r.sim.sediment: main program for hydrologic and sediment transport
* simulation (SIMWE)
*
* Original program (2002) and various modifications:
@@ -46,12 +46,19 @@
*
*/
+/********************************/
+/* DEFINE GLOB VAR */
+/********************************/
#define NWALK "2000000"
#define DIFFC "0.8"
#define NITER "1200"
#define ITEROUT "300"
#define DENSITY "200"
+#define MANINVAL "0.1"
+/********************************/
+/* INCLUDES */
+/********************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
@@ -67,6 +74,9 @@
#include <grass/site.h>
#include <grass/glocale.h>
+/********************************/
+/* Specific stuff */
+/********************************/
#define MAIN
#include <grass/waterglobs.h>
@@ -80,366 +90,344 @@
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, sediment flow, erosion, deposition");
+ module->description =
+ _("Sediment transport and erosion/deposition simulation "
+ "using path sampling method (SIMWE)");
- G_gisinit (argv[0]);
+ if (G_get_set_window (&cellhd) == -1)
+ exit (EXIT_FAILURE);
- 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();
- conv = G_database_units_to_meters_factor();
+ mixx = cellhd.west * conv;
+ maxx = cellhd.east * conv;
+ miyy = cellhd.south * conv;
+ mayy = cellhd.north * conv;
+
+ 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;
+
+ hhc = hhmax = 0.;
+
+ 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));
- mixx = cellhd.west * conv;
- maxx = cellhd.east * conv;
- miyy = cellhd.south * conv;
- mayy = cellhd.north * conv;
-
- 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;
+ 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");
- bxmi=2093113. * conv;
- bymi=731331. * conv;
- bxma=2093461. * conv;
- byma=731529. * conv;
- bresx=2. * conv;
- bresy=2. * conv;
- maxwab=100000;
+ parm.wdepth = G_define_standard_option(G_OPT_R_INPUT);
+ parm.wdepth->key = "wdepth";
+ parm.wdepth->description = _("Name of the water depth raster map [m]");
+ parm.wdepth->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.detin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.detin->key = "detin";
+ parm.detin->description =
+ _("Name of the detachment capacity coefficient raster map [s/m]");
+ parm.detin->guisection = _("Input_options");
+
+ parm.tranin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.tranin->key = "tranin";
+ parm.tranin->description =
+ _("Name of the transport capacity coefficient raster map [s]");
+ parm.tranin->guisection = _("Input_options");
+
+ parm.tauin = G_define_standard_option(G_OPT_R_INPUT);
+ parm.tauin->key = "tauin";
+ parm.tauin->description = _("Name of the critical shear stress raster map [Pa]");
+ parm.tauin->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 = "maninval";
+ parm.maninval->type = TYPE_DOUBLE;
+ parm.maninval->answer = MANINVAL;
+ parm.maninval->required = NO;
+ parm.maninval->description = _("Name of the Mannings n value");
+ parm.maninval->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.tc = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.tc->key = "tc";
+ parm.tc->required = NO;
+ parm.tc->description = _("Output transport capacity raster map [kg/ms]");
+ parm.tc->guisection = _("Output_options");
+
+ parm.et = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.et->key = "et";
+ parm.et->required = NO;
+ parm.et->description =
+ _("Output transp.limited erosion-deposition raster map [kg/m2s]");
+ parm.et->guisection = _("Output_options");
+
+ parm.conc = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.conc->key = "conc";
+ parm.conc->required = NO;
+ parm.conc->description = _("Output sediment concentration raster map [particle/m3]");
+ parm.conc->guisection = _("Output_options");
+
+ parm.flux = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.flux->key = "flux";
+ parm.flux->required = NO;
+ parm.flux->description = _("Output sediment flux raster map [kg/ms]");
+ parm.flux->guisection = _("Output_options");
+
+ parm.erdep = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.erdep->key = "erdep";
+ parm.erdep->required = NO;
+ parm.erdep->description = _("Output erosion-deposition raster map [kg/m2s]");
+ parm.erdep->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");
+ 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");
+
+ 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");
- hhc = hhmax = 0.;
+ if (G_parser (argc, argv))
+ exit (EXIT_FAILURE);
- mx2o= (int)((bxma-bxmi)/bresx);
- my2o= (int)((byma-bymi)/bresy);
+ mscale=flag.mscale->answer ? "m" : NULL;
+ tserie=flag.tserie->answer ? "t" : NULL;
+
+ elevin = parm.elevin->answer;
+ wdepth = parm.wdepth->answer;
+ dxin = parm.dxin->answer;
+ dyin = parm.dyin->answer;
+ /* maskmap = parm.maskmap->answer;*/
+ detin = parm.detin->answer;
+ tranin = parm.tranin->answer;
+ tauin = parm.tauin->answer;
+ manin = parm.manin->answer;
+ tc = parm.tc->answer;
+ et = parm.et->answer;
+ conc = parm.conc->answer;
+ flux = parm.flux->answer;
+ erdep = parm.erdep->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.maninval->answer, "%lf", &manin_val);
-/* relative small box coordinates: leave 1 grid layer for overlap */
+ /* 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 !!!!!"));
- 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));
+ /* 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;
+ }
+ /*rwalk = (double) maxwa;*/
- 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");
- parm.elevin->guisection = _("Input_options");
+ if (conv != 1.0)
+ G_message(_("Using metric conversion factor %f, step=%f"), conv, step);
- parm.wdepth = G_define_option();
- parm.wdepth->key = "wdepth";
- parm.wdepth->type = TYPE_STRING;
- parm.wdepth->required = YES;
- parm.wdepth->gisprompt = "old,cell,raster";
- parm.wdepth->description = _("Name of the water height raster map");
- parm.wdepth->guisection = _("Input_options");
+ /*
+ * G_set_embedded_null_value_mode(1);
+ */
+
+ if ((tc == NULL) && (et == NULL) && (conc == NULL) && (flux == NULL) && (erdep == NULL))
+ G_warning(_("You are not outputing any raster or site files"));
+ ret_val = input_data();
+ if (ret_val != 1)
+ G_fatal_error(_("Input failed"));
- 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");
- parm.dxin->guisection = _("Input_options");
+ /* mandatory for si,sigma */
- 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");
- parm.dyin->guisection = _("Input_options");
+ si = (double **)G_malloc (sizeof(double *)*(my));
+ for(l=0;l<my;l++)
+ {
+ si[l] = (double*)G_malloc (sizeof(double)*(mx));
+ }
+ for (j = 0; j < my; j++)
+ {
+ for (i = 0; i < mx; i++)
+ si[j][i] = 0.;
+ }
- parm.detin = G_define_option();
- parm.detin->key = "detin";
- parm.detin->type = TYPE_STRING;
- parm.detin->required = YES;
- parm.detin->gisprompt = "old,cell,raster";
- parm.detin->description =
- _("Name of the detachment capacity coefficient raster map");
- parm.detin->guisection = _("Input_options");
+ sigma = (double **)G_malloc (sizeof(double *)*(my));
+ for(l=0;l<my;l++)
+ {
+ sigma[l] = (double*)G_malloc (sizeof(double)*(mx));
+ }
+ for (j = 0; j < my; j++)
+ {
+ for (i = 0; i < mx; i++)
+ sigma[j][i] = 0.;
+ }
- parm.tranin = G_define_option();
- parm.tranin->key = "tranin";
- parm.tranin->type = TYPE_STRING;
- parm.tranin->required = YES;
- parm.tranin->gisprompt = "old,cell,raster";
- parm.tranin->description =
- _("Name of the transport capacity coefficient raster map");
- parm.tranin->guisection = _("Input_options");
+ /* memory allocation for output grids */
- parm.tauin = G_define_option();
- parm.tauin->key = "tauin";
- parm.tauin->type = TYPE_STRING;
- parm.tauin->required = YES;
- parm.tauin->gisprompt = "old,cell,raster";
- parm.tauin->description = _("Name of the critical shear stress raster map");
- parm.tauin->guisection = _("Input_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.;
+ }
- 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");
+ if (erdep != NULL || et != NULL)
+ {
+ er = (float **)G_malloc (sizeof(float *)*(my));
+ for(l=0;l<my;l++)
+ {
+ er[l] = (float*)G_malloc (sizeof(float)*(mx));
+ }
+ for (j = 0; j < my; j++)
+ {
+ for (i = 0; i < mx; i++)
+ er[j][i] = 0.;
+ }
+ }
-/* 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");
+ /* if (maskmap != NULL)
+ bitmask = BM_create (cols, rows);
+ IL_create_bitmask (¶ms, bitmask);
+ */
+ seeds(rand1,rand2);
+ grad_check();
- parm.tc = G_define_option();
- parm.tc->key = "tc";
- parm.tc->type = TYPE_STRING;
- parm.tc->required = NO;
- parm.tc->gisprompt = "new,cell,raster";
- parm.tc->description = _("Output transport capacity raster map");
- parm.tc->guisection = _("Output_options");
+ if (et != NULL)
+ erod(si);
+ /* treba dat output pre topoerdep */
+ main_loop();
- parm.et = G_define_option();
- parm.et->key = "et";
- parm.et->type = TYPE_STRING;
- parm.et->required = NO;
- parm.et->gisprompt = "new,cell,raster";
- parm.et->description =
- _("Output transp.limited erosion-deposition raster map");
- parm.et->guisection = _("Output_options");
+ if (tserie == NULL) {
+ ii=output_data (0,1.);
+ if (ii != 1)
+ G_fatal_error (_("Cannot write raster maps"));
+ }
- parm.conc = G_define_option();
- parm.conc->key = "conc";
- parm.conc->type = TYPE_STRING;
- parm.conc->required = NO;
- parm.conc->gisprompt = "new,cell,raster";
- parm.conc->description = _("Output sediment concentration raster map");
- parm.conc->guisection = _("Output_options");
+ if (fdwalkers != NULL)
+ fclose (fdwalkers);
- parm.flux = G_define_option();
- parm.flux->key = "flux";
- parm.flux->type = TYPE_STRING;
- parm.flux->required = NO;
- parm.flux->gisprompt = "new,cell,raster";
- parm.flux->description = _("Output sediment flux raster map");
- parm.flux->guisection = _("Output_options");
-
- parm.erdep = G_define_option();
- parm.erdep->key = "erdep";
- parm.erdep->type = TYPE_STRING;
- parm.erdep->required = NO;
- parm.erdep->gisprompt = "new,cell,raster";
- parm.erdep->description = _("Output erosion-deposition raster map");
- parm.erdep->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");
- 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");
-
- 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 ? "m" : NULL;
- tserie=flag.tserie->answer ? "t" : NULL;
-
- elevin = parm.elevin->answer;
- wdepth = parm.wdepth->answer;
- dxin = parm.dxin->answer;
- dyin = parm.dyin->answer;
-/* maskmap = parm.maskmap->answer;*/
- detin = parm.detin->answer;
- tranin = parm.tranin->answer;
- tauin = parm.tauin->answer;
- manin = parm.manin->answer;
- tc = parm.tc->answer;
- et = parm.et->answer;
- conc = parm.conc->answer;
- flux = parm.flux->answer;
- erdep = parm.erdep->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);
-
- 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 ((tc == NULL) && (et == NULL) && (conc == NULL) && (flux == NULL) && (erdep == NULL))
- G_warning(_("You are not outputing any raster or site files"));
-
- ret_val = input_data();
- if (ret_val != 1)
- G_fatal_error(_("Input failed"));
-
-/* mandatory for si,sigma */
-
-
- si = (double **)G_malloc (sizeof(double *)*(my));
- for(l=0;l<my;l++)
- {
- si[l] = (double*)G_malloc (sizeof(double)*(mx));
- }
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- si[j][i] = 0.;
- }
-
- sigma = (double **)G_malloc (sizeof(double *)*(my));
- for(l=0;l<my;l++)
- {
- sigma[l] = (double*)G_malloc (sizeof(double)*(mx));
- }
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- sigma[j][i] = 0.;
- }
-
-/* memory allocation for output grids */
-
- 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 (erdep != NULL || et != NULL)
- {
- er = (float **)G_malloc (sizeof(float *)*(my));
- for(l=0;l<my;l++)
- {
- er[l] = (float*)G_malloc (sizeof(float)*(mx));
- }
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- er[j][i] = 0.;
- }
- }
-
-
-
-
-
-/* if (maskmap != NULL)
- bitmask = BM_create (cols, rows);
-
- IL_create_bitmask (¶ms, bitmask);
-*/
-
-
- seeds(rand1,rand2);
-
- grad_check();
-
- if (et != NULL)
- erod(si);
-/* treba dat output pre topoerdep */
-
-/* fprintf (stdout, "\n Percent complete: ");*/
-
- main_loop();
-
- if (tserie == NULL) {
- ii=output_data (0,1.);
- if (ii != 1)
- G_fatal_error (_("Cannot write raster maps"));
- }
-
- if (fdwalkers != NULL)
- fclose (fdwalkers);
-
- if (sfile != NULL)
- G_sites_close(fw);
-
- exit (EXIT_SUCCESS);
+ if (sfile != NULL)
+ G_sites_close(fw);
+ /* Exit with Success */
+ exit (EXIT_SUCCESS);
}
Modified: grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/waterglobs.h
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/waterglobs.h 2008-02-22 12:44:58 UTC (rev 30285)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/waterglobs.h 2008-02-22 14:27:39 UTC (rev 30286)
@@ -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/r.sim.water/description.html
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/description.html 2008-02-22 12:44:58 UTC (rev 30285)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/description.html 2008-02-22 14:27:39 UTC (rev 30286)
@@ -61,8 +61,8 @@
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>
+Files are saved with a suffix representing time since the start of simulation in seconds
+(e.g. wdepth.500, wdepth.1000).<br>
<P>
Overland flow is routed based on partial derivatives of elevation
field or other landscape features influencing water flow. Simulation
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-22 12:44:58 UTC (rev 30285)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/main.c 2008-02-22 14:27:39 UTC (rev 30286)
@@ -379,7 +379,7 @@
}
}
/* Report the final value of rain_val*/
- G_message(_("rain_val is set to: %f\n"),rain_val);
+ G_debug(3, "rain_val is set to: %f\n", rain_val);
/* if no Mannings map, then:*/
if(parm.manin->answer==NULL){
@@ -405,7 +405,7 @@
}
}
/* Report the final value of manin_val*/
- G_message(_("manin_val is set to: %f\n"),manin_val);
+ G_debug(3,"manin_val is set to: %f\n",manin_val);
/* if no infiltration map, then:*/
if(parm.infil->answer==NULL){
@@ -431,7 +431,7 @@
}
}
/* Report the final value of infil_val*/
- G_message(_("infil_val is set to: %f\n"),infil_val);
+ G_debug(3,"infil_val is set to: %f\n",infil_val);
/* Recompute timesec from user input in minutes
* to real timesec in seconds */
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c 2008-02-22 12:44:58 UTC (rev 30285)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c 2008-02-22 14:27:39 UTC (rev 30286)
@@ -115,9 +115,9 @@
zz = (float **)G_malloc (sizeof(float *)*(my));
v1 = (double **)G_malloc (sizeof(double *)*(my));
v2 = (double **)G_malloc (sizeof(double *)*(my));
-
- if(rain != NULL||rain_val >= 0.0)
+ if(rain != NULL||rain_val >= 0.0){
si = (double **)G_malloc (sizeof(double *)*(my));
+ }
if(infil != NULL||infil_val >= 0.0)
inf = (double **)G_malloc (sizeof(double *)*(my));
@@ -169,7 +169,7 @@
gama[l] = (double*)G_malloc (sizeof(double)*(mx));
}
- G_message (_("Running MAY 10 version, started modifications on 20080211"));
+ G_debug(3, "Running MAY 10 version, started modifications on 20080211");
/* Check if data available in mapsets
* if found, then open the files */
@@ -280,7 +280,7 @@
for (j=0; j<mx; j++)
{
row_rev = my - row - 1;
-
+ /*if elevation data exists store in zz[][]*/
if(!G_is_f_null_value(cell1+j))
zz[row_rev][j] = (float ) (conv * cell1[j]);
else
@@ -296,32 +296,38 @@
else
v2[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 */
+ zz[row_rev][j] = UNDEF;
/* 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?
- */
+ * zz[row_rev][j] == UNDEF) {
+ * v1[row_rev][j] == UNDEF;
+ * v2[row_rev][j] == UNDEF;
+ * zz[row_rev][j] == UNDEF;
+ * }
+ */ /*printout warning?*/
+
+ /* If Rain Exists, then load data */
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*/
+ 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;
}
-
+
+ /* Load infiltration map too if it exists*/
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*/
+ 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;
@@ -348,43 +354,43 @@
}
}
} 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;
- }
+ /* 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){
@@ -394,7 +400,7 @@
cchez[row_rev][j] = UNDEF;
zz[row_rev][j] = UNDEF;
}
- } else if(manin_val>0.0) { /* Added by Yann 20080213 */
+ } 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);
@@ -682,7 +688,6 @@
} /*DEFined area */
}
}
-
return 1;
}
More information about the grass-commit
mailing list