[GRASS-SVN] r33426 - 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
Sat Sep 13 04:32:00 EDT 2008
Author: neteler
Date: 2008-09-13 04:32:00 -0400 (Sat, 13 Sep 2008)
New Revision: 33426
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/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/random.c
grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h
Log:
updated code layout to new rules
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-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/main.c 2008-09-13 08:32:00 UTC (rev 33426)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.sim.sediment: main program for sediment transport
@@ -48,8 +49,9 @@
/********************************/
/* DEFINE GLOB VAR */
+
/********************************/
-/* #define NWALK "1000000" */
+/* #define NWALK "1000000" */
#define DIFFC "0.8"
#define NITER "10"
#define ITEROUT "2"
@@ -58,6 +60,7 @@
/********************************/
/* INCLUDES */
+
/********************************/
#include <stdio.h>
#include <stdlib.h>
@@ -76,6 +79,7 @@
/********************************/
/* Specific stuff */
+
/********************************/
#define MAIN
#include <grass/waterglobs.h>
@@ -83,8 +87,8 @@
char fncdsm[32];
char filnam[10];
-/*struct BM *bitmask;*/
-/*struct Cell_head cellhd;*/
+/*struct BM *bitmask; */
+/*struct Cell_head cellhd; */
struct GModule *module;
struct Map_info Map;
@@ -92,344 +96,343 @@
/****************************************/
/* MAIN */
+
/****************************************/
-int main ( int argc, char *argv[])
+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;
-
- G_gisinit (argv[0]);
-
- module = G_define_module();
- module->keywords = _("raster, sediment flow, erosion, deposition");
- module->description =
+ 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)");
- if (G_get_set_window (&cellhd) == -1)
- exit (EXIT_FAILURE);
+ 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;
- hhc = hhmax = 0.;
-
-/*
- 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));
-*/
+ mixx = cellhd.west * conv;
+ maxx = cellhd.east * conv;
+ miyy = cellhd.south * conv;
+ mayy = cellhd.north * conv;
- 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.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");
-*/
+ 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;
+ hhc = hhmax = 0.;
- if (G_parser (argc, argv))
- exit (EXIT_FAILURE);
+ /*
+ bxmi=2093113. * conv;
+ bymi=731331. * conv;
+ bxma=2093461. * conv;
+ byma=731529. * conv;
+ bresx=2. * conv;
+ bresy=2. * conv;
+ maxwab=100000;
-/* 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);
+ mx2o= (int)((bxma-bxmi)/bresx);
+ my2o= (int)((byma-bymi)/bresy);
- /* 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 !!!!!"));
+ /* relative small box coordinates: leave 1 grid layer for overlap
- /* 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;*/
+ 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 (conv != 1.0)
- G_message(_("Using metric conversion factor %f, step=%f"), conv, step);
+ 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");
- /*
- * 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.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");
- /* mandatory for si,sigma */
+ 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");
- 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.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");
- 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.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");
- /* memory allocation for output grids */
+ 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");
- dif = (float **)G_malloc (sizeof(float *)*(my));
- for(l=0;l<my;l++)
- {
- dif[l] = (float*)G_malloc (sizeof(float)*(mx));
+ 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");
+ */
+
+ 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);
+ sscanf(parm.maninval->answer, "%lf", &manin_val);
+
+ /* 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 !!!!!"));
+
+ /* 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; */
+
+ 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++)
- dif[j][i] = 0.;
+ for (j = 0; j < my; j++) {
+ for (i = 0; i < mx; i++)
+ er[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 (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 */
+ main_loop();
- if (et != NULL)
- erod(si);
- /* treba dat output pre topoerdep */
- main_loop();
+ if (tserie == NULL) {
+ ii = output_data(0, 1.);
+ if (ii != 1)
+ G_fatal_error(_("Cannot write raster maps"));
+ }
- if (tserie == NULL) {
- ii=output_data (0,1.);
- if (ii != 1)
- G_fatal_error (_("Cannot write raster maps"));
- }
+ if (fdwalkers != NULL)
+ fclose(fdwalkers);
- if (fdwalkers != NULL)
- fclose (fdwalkers);
-
- if (sfile != NULL)
- G_sites_close(fw);
- /* Exit with Success */
- 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-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.sediment/waterglobs.h 2008-09-13 08:32:00 UTC (rev 33426)
@@ -3,9 +3,11 @@
#define UNDEF -9999
#ifdef MAIN
-FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps, *fdmanin, *fddepth, *fddisch, *fderr, *fdoutwalk,*fdwalkers;
-FILE *fdwdepth,*fddetin,*fdtranin,*fdtauin, *fdtc, *fdet, *fdconc, *fdflux,*fderdep;
-FILE *fdsfile,*fw;
+FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps, *fdmanin,
+ *fddepth, *fddisch, *fderr, *fdoutwalk, *fdwalkers;
+FILE *fdwdepth, *fddetin, *fdtranin, *fdtauin, *fdtc, *fdet, *fdconc, *fdflux,
+ *fderdep;
+FILE *fdsfile, *fw;
char *elevin;
char *dxin;
@@ -21,7 +23,7 @@
char *outwalk;
char *mapset;
-/* Flag*/
+/* Flag */
char *tserie;
char *wdepth;
@@ -34,25 +36,28 @@
char *flux;
char *erdep;
-/*New Input val*/
+/*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,*rainval,*maninval,*infilval;
- } parm;
+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, *rainval, *maninval,
+ *infilval;
+} parm;
- struct
- {
+struct
+{
struct Flag *tserie;
- } flag;
+} flag;
-struct {
+struct
+{
long int is1, is2;
} seed;
@@ -69,53 +74,53 @@
int input_data(void);
-int seeds(long int,long int);
-int seedg(long int,long int);
+int seeds(long int, long int);
+int seedg(long int, long int);
int grad_check(void);
void erod(double **);
void main_loop(void);
-int output_data(int,double);
+int output_data(int, double);
int output_et(void);
double ulec(void);
double gasdev(void);
-double amax1(double,double);
-double amin1(double,double);
-int min(int,int);
-int max(int,int);
+double amax1(double, double);
+double amin1(double, double);
+int min(int, int);
+int max(int, int);
double xmin, ymin, xmax, ymax;
double mayy, miyy, maxx, mixx;
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;
+double step, conv;
double frac;
double bxmi, bymi;
float **zz, **cchez;
double **v1, **v2, **slope;
-double **gama, **gammas,**si,**inf,**sigma;
-float **dc,**tau,**er, **ct, **trap;
-float **dif;
+double **gama, **gammas, **si, **inf, **sigma;
+float **dc, **tau, **er, **ct, **trap;
+float **dif;
-double vavg[MAXW][2], stack[MAXW][3],w[MAXW][3];
+double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
int iflag[MAXW];
double hbeta;
int ldemo;
-double hhmax, sisum,vmean;
-double infsum,infmean;
+double hhmax, sisum, vmean;
+double infsum, infmean;
int maxw, maxwa, nwalk;
double rwalk, bresx, bresy, xrand, yrand;
double stepx, stepy, xp0, yp0;
-double chmean, si0, deltap, deldif, cch, hhc,halpha;
+double chmean, si0, deltap, deldif, cch, hhc, halpha;
double eps;
int maxwab, nstack;
int iterout, mx2o, my2o;
-int miter,nwalka,lwwfin;
+int miter, nwalka, lwwfin;
double timec;
int ts, timesec;
@@ -124,9 +129,11 @@
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;
-extern FILE *fdsfile,*fw;
+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;
+extern FILE *fdsfile, *fw;
extern char *elevin;
extern char *dxin;
@@ -158,19 +165,22 @@
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,*rainval,*maninval,*infilval;
- } parm;
+{
+ 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, *rainval, *maninval,
+ *infilval;
+} parm;
extern struct
- {
+{
struct Flag *tserie;
- } flag;
+} flag;
-extern struct {
+extern struct
+{
long int is1, is2;
} seed;
@@ -188,53 +198,53 @@
extern int input_data(void);
-extern int seeds(long int,long int);
-extern int seedg(long int,long int);
+extern int seeds(long int, long int);
+extern int seedg(long int, long int);
extern int grad_check(void);
extern void erod(double **);
extern void main_loop(void);
-extern int output_data(int,double);
+extern int output_data(int, double);
extern int output_et(void);
extern double ulec(void);
extern double gasdev(void);
-extern double amax1(double,double);
-extern double amin1(double,double);
-extern int min(int,int);
-extern int max(int,int);
+extern double amax1(double, double);
+extern double amin1(double, double);
+extern int min(int, int);
+extern int max(int, int);
extern double xmin, ymin, xmax, ymax;
extern double mayy, miyy, maxx, mixx;
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;
+extern double step, conv;
extern double frac;
extern double bxmi, bymi;
extern float **zz, **cchez;
extern double **v1, **v2, **slope;
-extern double **gama, **gammas,**si,**inf,**sigma;
-extern float **dc,**tau,**er, **ct, **trap;
-extern float **dif;
+extern double **gama, **gammas, **si, **inf, **sigma;
+extern float **dc, **tau, **er, **ct, **trap;
+extern float **dif;
-extern double vavg[MAXW][2], stack[MAXW][3],w[MAXW][3];
+extern double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
extern int iflag[MAXW];
extern double hbeta;
extern int ldemo;
-extern double hhmax, sisum,vmean;
-extern double infsum,infmean;
+extern double hhmax, sisum, vmean;
+extern double infsum, infmean;
extern int maxw, maxwa, nwalk;
extern double rwalk, bresx, bresy, xrand, yrand;
extern double stepx, stepy, xp0, yp0;
-extern double chmean, si0, deltap, deldif, cch, hhc,halpha;
+extern double chmean, si0, deltap, deldif, cch, hhc, halpha;
extern double eps;
extern int maxwab, nstack;
extern int iterout, mx2o, my2o;
-extern int miter,nwalka,lwwfin;
+extern int miter, nwalka, lwwfin;
extern double timec;
extern int ts, timesec;
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-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/main.c 2008-09-13 08:32:00 UTC (rev 33426)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.sim.water: main program for hydrologic and sediment transport
@@ -47,9 +48,10 @@
*/
/********************************/
-/* DEFINE GLOB VAR */
+/* DEFINE GLOB VAR */
+
/********************************/
-/* #define NWALK "1000000" */
+/* #define NWALK "1000000" */
#define DIFFC "0.8"
#define HMAX "0.3"
#define HALPHA "4.0"
@@ -62,7 +64,8 @@
#define INFILVAL "0.0"
/********************************/
-/* INCLUDES */
+/* INCLUDES */
+
/********************************/
#include <stdio.h>
#include <stdlib.h>
@@ -80,7 +83,8 @@
#include <grass/glocale.h>
/********************************/
-/* Specific stuff */
+/* Specific stuff */
+
/********************************/
#define MAIN
#include <grass/waterglobs.h>
@@ -88,441 +92,457 @@
char fncdsm[32];
char filnam[10];
-/*struct BM *bitmask;*/
-/*struct Cell_head cellhd;*/
+/*struct BM *bitmask; */
+/*struct Cell_head cellhd; */
struct GModule *module;
struct Map_info Map;
char msg[1024];
/****************************************/
-/* MAIN */
+/* MAIN */
+
/****************************************/
-int main ( int argc, char *argv[])
+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;
-
- 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);
+ int i, ii, j, l;
+ int ret_val;
+ double x_orig, y_orig;
+ static int rand1 = 12345;
+ static int rand2 = 67891;
- 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");
+ 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_parser (argc, argv))
- exit (EXIT_FAILURE);
+ if (G_get_set_window(&cellhd) == -1)
+ 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);
+ conv = G_database_units_to_meters_factor();
- /* 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;
- }
+ 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");
+
+
+ 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);
+
+ /* 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 */
}
- /* Report the final value of rain_val*/
- G_debug(3, "rain_val is set to: %f\n", rain_val);
+ 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_debug(3, "rain_val is set to: %f\n", rain_val);
- /* 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;
- }
+ /* 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 */
}
- /* Report the final value of manin_val*/
- G_debug(3,"manin_val is set to: %f\n",manin_val);
+ 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_debug(3, "manin_val is set to: %f\n", manin_val);
- /* 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 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;
- }
+ }
+ else {
+ /* Map input, set infil_val to -999.99 */
+ if (parm.infilval->answer == NULL) {
+ infil_val = -999.99;
}
- /* Report the final value of infil_val*/
- G_debug(3,"infil_val is set to: %f\n",infil_val);
+ 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_debug(3, "infil_val is set to: %f\n", infil_val);
- /* 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 !!!!!"));
+ /* 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 !!!!!"));
- /* 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;
- }
+ /* 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; */
+ /* rwalk = (double) maxwa; */
- if (conv != 1.0)
- G_message(_("Using metric conversion factor %f, step=%f"), conv, step);
+ if (conv != 1.0)
+ G_message(_("Using metric conversion factor %f, step=%f"), conv,
+ step);
- /*
- * G_set_embedded_null_value_mode(1);
- */
+ /*
+ * 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"));
+ 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 */
+ /* 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.;
- }
+ 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.;
- }
+ if (err != NULL) {
+ gammas = (double **)G_malloc(sizeof(double *) * (my));
+ for (l = 0; l < my; l++) {
+ gammas[l] = (double *)G_malloc(sizeof(double) * (mx));
}
-
- 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++)
+ gammas[j][i] = 0.;
}
- 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();
+ }
- if (ts == 0) {
- ii=output_data (0,1.);
- if (ii != 1)
- G_fatal_error(_("Cannot write raster maps"));
- }
+ 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 (fdwalkers != NULL)
- fclose (fdwalkers);
+ /* if (maskmap != NULL)
+ bitmask = BM_create (cols, rows);
+ IL_create_bitmask (¶ms, bitmask);
+ */
+ seeds(rand1, rand2);
+ grad_check();
+ main_loop();
- if (sfile != NULL)
- fclose(fw);
- /* Exit with Success */
- exit (EXIT_SUCCESS);
+ 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 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-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/r.sim.water/waterglobs.h 2008-09-13 08:32:00 UTC (rev 33426)
@@ -3,9 +3,11 @@
#define UNDEF -9999
#ifdef MAIN
-FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps, *fdmanin, *fddepth, *fddisch, *fderr, *fdoutwalk,*fdwalkers;
-FILE *fdwdepth,*fddetin,*fdtranin,*fdtauin, *fdtc, *fdet, *fdconc, *fdflux,*fderdep;
-FILE *fdsfile,*fw;
+FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps, *fdmanin,
+ *fddepth, *fddisch, *fderr, *fdoutwalk, *fdwalkers;
+FILE *fdwdepth, *fddetin, *fdtranin, *fdtauin, *fdtc, *fdet, *fdconc, *fdflux,
+ *fderdep;
+FILE *fdsfile, *fw;
char *elevin;
char *dxin;
@@ -21,7 +23,7 @@
char *outwalk;
char *mapset;
-/* Flag*/
+/* Flag */
char *tserie;
char *wdepth;
@@ -34,25 +36,28 @@
char *flux;
char *erdep;
-/*New Input val*/
+/*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,*rainval,*maninval,*infilval;
- } parm;
+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, *rainval, *maninval,
+ *infilval;
+} parm;
- struct
- {
+struct
+{
struct Flag *tserie;
- } flag;
+} flag;
-struct {
+struct
+{
long int is1, is2;
} seed;
@@ -69,53 +74,53 @@
int input_data(void);
-int seeds(long int,long int);
-int seedg(long int,long int);
+int seeds(long int, long int);
+int seedg(long int, long int);
int grad_check(void);
void erod(double **);
void main_loop(void);
-int output_data(int,double);
+int output_data(int, double);
int output_et(void);
double ulec(void);
double gasdev(void);
-double amax1(double,double);
-double amin1(double,double);
-int min(int,int);
-int max(int,int);
+double amax1(double, double);
+double amin1(double, double);
+int min(int, int);
+int max(int, int);
double xmin, ymin, xmax, ymax;
double mayy, miyy, maxx, mixx;
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;
+double step, conv;
double frac;
double bxmi, bymi;
float **zz, **cchez;
double **v1, **v2, **slope;
-double **gama, **gammas,**si,**inf,**sigma;
-float **dc,**tau,**er, **ct, **trap;
-float **dif;
+double **gama, **gammas, **si, **inf, **sigma;
+float **dc, **tau, **er, **ct, **trap;
+float **dif;
-double vavg[MAXW][2], stack[MAXW][3],w[MAXW][3];
+double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
int iflag[MAXW];
double hbeta;
int ldemo;
-double hhmax, sisum,vmean;
-double infsum,infmean;
+double hhmax, sisum, vmean;
+double infsum, infmean;
int maxw, maxwa, nwalk;
double rwalk, bresx, bresy, xrand, yrand;
double stepx, stepy, xp0, yp0;
-double chmean, si0, deltap, deldif, cch, hhc,halpha;
+double chmean, si0, deltap, deldif, cch, hhc, halpha;
double eps;
int maxwab, nstack;
int iterout, mx2o, my2o;
-int miter,nwalka,lwwfin;
+int miter, nwalka, lwwfin;
double timec;
int ts, timesec;
@@ -124,9 +129,11 @@
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;
-extern FILE *fdsfile,*fw;
+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;
+extern FILE *fdsfile, *fw;
extern char *elevin;
extern char *dxin;
@@ -158,19 +165,22 @@
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,*rainval,*maninval,*infilval;
- } parm;
+{
+ 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, *rainval, *maninval,
+ *infilval;
+} parm;
extern struct
- {
+{
struct Flag *tserie;
- } flag;
+} flag;
-extern struct {
+extern struct
+{
long int is1, is2;
} seed;
@@ -188,53 +198,53 @@
extern int input_data(void);
-extern int seeds(long int,long int);
-extern int seedg(long int,long int);
+extern int seeds(long int, long int);
+extern int seedg(long int, long int);
extern int grad_check(void);
extern void erod(double **);
extern void main_loop(void);
-extern int output_data(int,double);
+extern int output_data(int, double);
extern int output_et(void);
extern double ulec(void);
extern double gasdev(void);
-extern double amax1(double,double);
-extern double amin1(double,double);
-extern int min(int,int);
-extern int max(int,int);
+extern double amax1(double, double);
+extern double amin1(double, double);
+extern int min(int, int);
+extern int max(int, int);
extern double xmin, ymin, xmax, ymax;
extern double mayy, miyy, maxx, mixx;
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;
+extern double step, conv;
extern double frac;
extern double bxmi, bymi;
extern float **zz, **cchez;
extern double **v1, **v2, **slope;
-extern double **gama, **gammas,**si,**inf,**sigma;
-extern float **dc,**tau,**er, **ct, **trap;
-extern float **dif;
+extern double **gama, **gammas, **si, **inf, **sigma;
+extern float **dc, **tau, **er, **ct, **trap;
+extern float **dif;
-extern double vavg[MAXW][2], stack[MAXW][3],w[MAXW][3];
+extern double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
extern int iflag[MAXW];
extern double hbeta;
extern int ldemo;
-extern double hhmax, sisum,vmean;
-extern double infsum,infmean;
+extern double hhmax, sisum, vmean;
+extern double infsum, infmean;
extern int maxw, maxwa, nwalk;
extern double rwalk, bresx, bresy, xrand, yrand;
extern double stepx, stepy, xp0, yp0;
-extern double chmean, si0, deltap, deldif, cch, hhc,halpha;
+extern double chmean, si0, deltap, deldif, cch, hhc, halpha;
extern double eps;
extern int maxwab, nstack;
extern int iterout, mx2o, my2o;
-extern int miter,nwalka,lwwfin;
+extern int miter, nwalka, lwwfin;
extern double timec;
extern int ts, timesec;
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/erod.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/erod.c 2008-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/erod.c 2008-09-13 08:32:00 UTC (rev 33426)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: simwe library
@@ -28,41 +29,41 @@
void erod(double **hw)
{
-/* hw = sigma or gamma */
+ /* hw = sigma or gamma */
- double dyp, dyn, dya, dxp, dxn, dxa;
- int k,l;
- int l1, lp, k1, kp, ln, kn, k2, l2;
+ double dyp, dyn, dya, dxp, dxn, dxa;
+ int k, l;
+ int l1, lp, k1, kp, ln, kn, k2, l2;
-
- for (k = 0; k < my; k++) {
- for (l = 0; l < mx; l++) {
- lp = max(0,l-2);
- l1 = lp + 1;
- kp = max(0,k-2);
- k1 = kp + 1;
- ln = min(mx-1,l+1);
- l2 = ln - 1;
- kn = min(my-1,k+1);
- k2 = kn - 1;
- if (zz[k][l] != UNDEF || zz[k][ln] != UNDEF || zz[kp][l] != UNDEF || \
- zz[k][lp] != UNDEF || zz[k][l1] != UNDEF || \
- zz[k1][l] != UNDEF || zz[kn][l] != UNDEF) { /* jh fix */
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
+ lp = max(0, l - 2);
+ l1 = lp + 1;
+ kp = max(0, k - 2);
+ k1 = kp + 1;
+ ln = min(mx - 1, l + 1);
+ l2 = ln - 1;
+ kn = min(my - 1, k + 1);
+ k2 = kn - 1;
+ if (zz[k][l] != UNDEF || zz[k][ln] != UNDEF || zz[kp][l] != UNDEF || zz[k][lp] != UNDEF || zz[k][l1] != UNDEF || zz[k1][l] != UNDEF || zz[kn][l] != UNDEF) { /* jh fix */
+
dxp = (v1[k][lp] * hw[k][lp] - v1[k][l1] * hw[k][l1]) / stepx;
- dxn = (v1[k][l2] * hw[k][l2] - v1[k][ln] * hw[k][ln]) / stepx;
+ dxn = (v1[k][l2] * hw[k][l2] - v1[k][ln] * hw[k][ln]) / stepx;
dxa = 0.5 * (dxp + dxn);
- dyp = (v2[kp][l] * hw[kp][l] - v2[k1][l] * hw[k1][l]) / stepy;
- dyn = (v2[k2][l] * hw[k2][l] - v2[kn][l] * hw[kn][l]) / stepy;
- dya = 0.5 * (dyp + dyn);
+ dyp = (v2[kp][l] * hw[kp][l] - v2[k1][l] * hw[k1][l]) / stepy;
+ dyn = (v2[k2][l] * hw[k2][l] - v2[kn][l] * hw[kn][l]) / stepy;
+ dya = 0.5 * (dyp + dyn);
er[k][l] = (dxa + dya) / deltap;
- } else er[k][l] = UNDEF;
+ }
+ else
+ er[k][l] = UNDEF;
- }
}
+ }
}
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/hydro.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/hydro.c 2008-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/hydro.c 2008-09-13 08:32:00 UTC (rev 33426)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: simwe library
@@ -31,38 +32,38 @@
/* ******************************************************** */
/* .......... iblock loop */
-void main_loop (void)
+void main_loop(void)
{
- int i,ii,l,k;
- int icoub,lwout,nmult;
+ int i, ii, l, k;
+ int icoub, lwout, nmult;
int iw, iblock, lw;
int itime, iter1;
- int nfiterh,nfiterw;
+ int nfiterh, nfiterw;
int mgen, mgen2, mgen3;
int nblock;
int icfl;
- int mitfac,p;
+ int mitfac, p;
double x, y;
double velx, vely, stxm, stym;
- double factor,conn,gaux, gauy;
- double d1,addac, decr;
- double barea,sarea,walkwe;
- double gen, gen2,wei2, wei3,wei,weifac;
+ double factor, conn, gaux, gauy;
+ double d1, addac, decr;
+ double barea, sarea, walkwe;
+ double gen, gen2, wei2, wei3, wei, weifac;
float eff;
nblock = 1;
icoub = 0;
icfl = 0;
nstack = 0;
-
- if (maxwa > (MAXW - mx * my)) {
+
+ if (maxwa > (MAXW - mx * my)) {
mitfac = maxwa / (MAXW - mx * my);
nblock = mitfac + 1;
maxwa = maxwa / nblock;
- }
+ }
-G_debug(2, " maxwa, nblock %d %d", maxwa, nblock);
+ G_debug(2, " maxwa, nblock %d %d", maxwa, nblock);
for (iblock = 1; iblock <= nblock; iblock++) {
++icoub;
@@ -71,313 +72,324 @@
walkwe = 0.;
barea = stepx * stepy;
sarea = bresx * bresy;
-G_debug(2, " barea,sarea,rwalk,sisum: %f %f %f %f",barea,sarea,rwalk,sisum);
+ G_debug(2, " barea,sarea,rwalk,sisum: %f %f %f %f", barea, sarea,
+ rwalk, sisum);
lwwfin = 0;
-/* write hh.walkers0 */
+ /* write hh.walkers0 */
for (k = 0; k < my; k++) {
- for (l = 0; l < mx; l++) { /* run thru the whole area */
- if (zz[k][l] != UNDEF){
+ for (l = 0; l < mx; l++) { /* run thru the whole area */
+ if (zz[k][l] != UNDEF) {
- x = xp0 + stepx * (double) (l);
- y = yp0 + stepy * (double) (k);
+ x = xp0 + stepx * (double)(l);
+ y = yp0 + stepy * (double)(k);
- gen = rwalk * si[k][l] / sisum;
- mgen = (int) gen;
- wei = gen / (double) (mgen + 1);
+ gen = rwalk * si[k][l] / sisum;
+ mgen = (int)gen;
+ wei = gen / (double)(mgen + 1);
- /*if (si[k][l] != 0.) {*/
-/* this stuff later for multiscale */
+ /*if (si[k][l] != 0.) { */
+ /* this stuff later for multiscale */
- gen2=(double)maxwab*si[k][l]/(si0*(double)(mx2o * my2o));
+ gen2 =
+ (double)maxwab *si[k][l] / (si0 *
+ (double)(mx2o * my2o));
gen2 = gen2 * (barea / sarea);
- mgen2 = (int) gen2;
- wei2 = gen2 / (double) (mgen2 + 1);
- mgen3 = (int) ((double) mgen2 * wei2 / ((double) mgen * wei));
+ mgen2 = (int)gen2;
+ wei2 = gen2 / (double)(mgen2 + 1);
+ mgen3 =
+ (int)((double)mgen2 * wei2 / ((double)mgen * wei));
nmult = mgen3 + 1;
- wei3 = gen2 / (double) ((mgen + 1) * (mgen2 + 1));
+ wei3 = gen2 / (double)((mgen + 1) * (mgen2 + 1));
weifac = wei3 / wei;
-/* } else {
- nmult = 1;
- weifac = 1.;
- fprintf(stderr, "\n zero rainfall excess in cell");
- }*/
+ /* } else {
+ nmult = 1;
+ weifac = 1.;
+ fprintf(stderr, "\n zero rainfall excess in cell");
+ } */
-/*G_debug(2, " gen,gen2,wei,wei2,mgen3,nmult: %f %f %f %f %d %d",gen,gen2,wei,wei2,mgen3,nmult);
-*/
- for (iw = 1; iw <= mgen+1; iw++) { /* assign walkers */
- ++lw;
+ /*G_debug(2, " gen,gen2,wei,wei2,mgen3,nmult: %f %f %f %f %d %d",gen,gen2,wei,wei2,mgen3,nmult);
+ */
+ for (iw = 1; iw <= mgen + 1; iw++) { /* assign walkers */
+ ++lw;
- if (lw > MAXW)
- G_fatal_error (_("nwalk > maxw!"));
+ if (lw > MAXW)
+ G_fatal_error(_("nwalk > maxw!"));
- w[lw][1] = x + stepx * (ulec() - 0.5);
- w[lw][2] = y + stepy * (ulec() - 0.5);
- w[lw][3] = wei;
+ w[lw][1] = x + stepx * (ulec() - 0.5);
+ w[lw][2] = y + stepy * (ulec() - 0.5);
+ w[lw][3] = wei;
- walkwe += w[lw][3];
- vavg[lw][1] = v1[k][l];
- vavg[lw][2] = v2[k][l];
- if(w[lw][1] >= xmin && w[lw][2] >= ymin && w[lw][1] <= xmax && w[lw][2] <= ymax)
- {
- iflag[lw] = 0;
- } else {
- iflag[lw] = 1;
- }
+ walkwe += w[lw][3];
+ vavg[lw][1] = v1[k][l];
+ vavg[lw][2] = v2[k][l];
+ if (w[lw][1] >= xmin && w[lw][2] >= ymin &&
+ w[lw][1] <= xmax && w[lw][2] <= ymax) {
+ iflag[lw] = 0;
+ }
+ else {
+ iflag[lw] = 1;
+ }
- lwout = lw / ldemo;
- lwout *= ldemo;
+ lwout = lw / ldemo;
+ lwout *= ldemo;
- if (lwout == lw) {
- ++lwwfin;
+ if (lwout == lw) {
+ ++lwwfin;
+ }
}
- }
- } /*DEFined area */
+ } /*DEFined area */
}
}
nwalk = lw;
-G_debug(2, " number of written walkers: %d",lwwfin);
-G_debug(2, " nwalk, maxw %d %d", nwalk, MAXW);
-G_debug(2, " walkwe (walk weight),frac %f %f",walkwe,frac);
+ G_debug(2, " number of written walkers: %d", lwwfin);
+ G_debug(2, " nwalk, maxw %d %d", nwalk, MAXW);
+ G_debug(2, " walkwe (walk weight),frac %f %f", walkwe, frac);
- stxm = stepx * (double) (mx+1) - xmin;
- stym = stepy * (double) (my+1) - ymin;
+ stxm = stepx * (double)(mx + 1) - xmin;
+ stym = stepy * (double)(my + 1) - ymin;
nwalka = 0;
- deldif = sqrt(deltap) * frac; /* diffuse factor */
+ deldif = sqrt(deltap) * frac; /* diffuse factor */
- factor = deltap*sisum/(rwalk * (double) nblock);
+ factor = deltap * sisum / (rwalk * (double)nblock);
-G_debug(2, " deldif,factor %f %e",deldif,factor);
+ G_debug(2, " deldif,factor %f %e", deldif, factor);
-/* ********************************************************** */
-/* main loop over the projection time */
-/* *********************************************************** */
+ /* ********************************************************** */
+ /* main loop over the projection time */
+ /* *********************************************************** */
-G_debug(2, "main loop over the projection time... ");
+ G_debug(2, "main loop over the projection time... ");
- for (i = 1; i <= miter; i++) { /* iteration loop depending on simulation time and deltap */
- G_percent(i,miter,1);
+ for (i = 1; i <= miter; i++) { /* iteration loop depending on simulation time and deltap */
+ G_percent(i, miter, 1);
iter1 = i / iterout;
iter1 *= iterout;
if (iter1 == i) {
nfiterw = i / iterout + 10;
nfiterh = i / iterout + 40;
- G_debug(2, "iblock=%d i=%d miter=%d nwalk=%d nwalka=%d",iblock,i,miter,nwalk,nwalka);
+ G_debug(2, "iblock=%d i=%d miter=%d nwalk=%d nwalka=%d",
+ iblock, i, miter, nwalk, nwalka);
}
-
- if (nwalka == 0 && i > 1) goto L_800;
-/* ************************************************************ */
-/* .... propagate one step */
-/* ************************************************************ */
+ if (nwalka == 0 && i > 1)
+ goto L_800;
+ /* ************************************************************ */
+ /* .... propagate one step */
+ /* ************************************************************ */
+
addac = factor;
- conn = (double) nblock / (double) iblock;
+ conn = (double)nblock / (double)iblock;
if (i == 1) {
addac = factor * .5;
}
nwalka = 0;
for (lw = 1; lw <= nwalk; lw++) {
- if (w[lw][3] > EPS) { /* check the walker weight */
+ if (w[lw][3] > EPS) { /* check the walker weight */
++nwalka;
- l = (int) ((w[lw][1] + stxm) / stepx) - mx -1;
- k = (int) ((w[lw][2] + stym) / stepy) - my -1;
+ l = (int)((w[lw][1] + stxm) / stepx) - mx - 1;
+ k = (int)((w[lw][2] + stym) / stepy) - my - 1;
- if (l > mx-1 || k > my-1 || k < 0 || l < 0) {
+ if (l > mx - 1 || k > my - 1 || k < 0 || l < 0) {
- G_debug(2, " k,l=%d,%d",k,l);
- printf(" lw,w=%d %f %f",lw,w[lw][1],w[lw][2]);
- G_debug(2, " stxym=%f %f",stxm,stym);
- printf(" step=%f %f",stepx,stepy);
- G_debug(2, " m=%d %d",my,mx);
- printf(" nwalka,nwalk=%d %d",nwalka,nwalk);
- G_debug(2, " ");
+ G_debug(2, " k,l=%d,%d", k, l);
+ printf(" lw,w=%d %f %f", lw, w[lw][1], w[lw][2]);
+ G_debug(2, " stxym=%f %f", stxm, stym);
+ printf(" step=%f %f", stepx, stepy);
+ G_debug(2, " m=%d %d", my, mx);
+ printf(" nwalka,nwalk=%d %d", nwalka, nwalk);
+ G_debug(2, " ");
}
if (zz[k][l] != UNDEF) {
- if(infil!=NULL) { /* infiltration part */
+ if (infil != NULL) { /* infiltration part */
- if (inf[k][l] - si[k][l] > 0.) {
+ if (inf[k][l] - si[k][l] > 0.) {
- decr = pow(addac * w[lw][3],3./5.); /* decreasing factor in m */
- if (inf[k][l] > decr) {
- inf[k][l] -= decr; /* decrease infilt. in cell and eliminate the walker */
- w[lw][3] = 0.;
- }
- else {
- w[lw][3] -= pow(inf[k][l],5./3.) / addac; /* use just proportional part of the walker weight */
- inf[k][l] = 0.;
+ decr = pow(addac * w[lw][3], 3. / 5.); /* decreasing factor in m */
+ if (inf[k][l] > decr) {
+ inf[k][l] -= decr; /* decrease infilt. in cell and eliminate the walker */
+ w[lw][3] = 0.;
+ }
+ else {
+ w[lw][3] -= pow(inf[k][l], 5. / 3.) / addac; /* use just proportional part of the walker weight */
+ inf[k][l] = 0.;
- }
+ }
- }
+ }
- }
+ }
- gama[k][l] += (addac * w[lw][3]); /* add walker weigh to water depth or conc. */
+ gama[k][l] += (addac * w[lw][3]); /* add walker weigh to water depth or conc. */
- d1 = gama[k][l] * conn;
- hhc = pow(d1, 3./5.);
+ d1 = gama[k][l] * conn;
+ hhc = pow(d1, 3. / 5.);
- if (hhc > hhmax && wdepth == NULL) { /* increased diffusion if w.depth > hhmax */
- dif[k][l] = (halpha + 1) * deldif;
- velx = vavg[lw][1];
- vely = vavg[lw][2];
- } else {
- dif[k][l] = deldif;
- velx = v1[k][l];
- vely = v2[k][l];
- }
+ if (hhc > hhmax && wdepth == NULL) { /* increased diffusion if w.depth > hhmax */
+ dif[k][l] = (halpha + 1) * deldif;
+ velx = vavg[lw][1];
+ vely = vavg[lw][2];
+ }
+ else {
+ dif[k][l] = deldif;
+ velx = v1[k][l];
+ vely = v2[k][l];
+ }
- if (traps != NULL && trap[k][l] != 0.) { /* traps */
+ if (traps != NULL && trap[k][l] != 0.) { /* traps */
- eff = ulec(); /* random generator */
+ eff = ulec(); /* random generator */
- if(eff <= trap[k][l]) {
- velx = - 0.1 * v1[k][l]; /* move it slightly back */
- vely = - 0.1 * v2[k][l];
+ if (eff <= trap[k][l]) {
+ velx = -0.1 * v1[k][l]; /* move it slightly back */
+ vely = -0.1 * v2[k][l];
+ }
+
}
- }
+ gaux = gasdev();
+ gauy = gasdev();
- gaux = gasdev();
- gauy = gasdev();
+ w[lw][1] += (velx + dif[k][l] * gaux); /* move the walker */
+ w[lw][2] += (vely + dif[k][l] * gauy);
- w[lw][1] += (velx + dif[k][l] * gaux); /* move the walker */
- w[lw][2] += (vely + dif[k][l] * gauy);
+ if (hhc > hhmax && wdepth == NULL) {
+ vavg[lw][1] = hbeta * (vavg[lw][1] + v1[k][l]);
+ vavg[lw][2] = hbeta * (vavg[lw][2] + v2[k][l]);
+ }
- if (hhc > hhmax && wdepth == NULL) {
- vavg[lw][1] = hbeta * (vavg[lw][1] + v1[k][l]);
- vavg[lw][2] = hbeta * (vavg[lw][2] + v2[k][l]);
- }
-
- if (w[lw][1] <= xmin || w[lw][2] <= ymin || w[lw][1]
+ if (w[lw][1] <= xmin || w[lw][2] <= ymin || w[lw][1]
>= xmax || w[lw][2] >= ymax) {
- w[lw][3] = 1e-10; /* eliminate walker if it is out of area */
- } else
- {
- if(wdepth!=NULL) {
- l = (int) ((w[lw][1] + stxm) / stepx) - mx -1;
- k = (int) ((w[lw][2] + stym) / stepy) - my -1;
+ w[lw][3] = 1e-10; /* eliminate walker if it is out of area */
+ }
+ else {
+ if (wdepth != NULL) {
+ l = (int)((w[lw][1] + stxm) / stepx) - mx - 1;
+ k = (int)((w[lw][2] + stym) / stepy) - my - 1;
w[lw][3] *= sigma[k][l];
- }
+ }
- } /* else*/
- } /*DEFined area */
+ } /* else */
+ } /*DEFined area */
else {
- w[lw][3] = 1e-10; /* eliminate walker if it is out of area */
+ w[lw][3] = 1e-10; /* eliminate walker if it is out of area */
}
}
-/* write the walkers which reach the box */
+ /* write the walkers which reach the box */
- if (w[lw][1] >= xmin && w[lw][2] >= ymin && w[lw][1]
- <= xmax && w[lw][2] <= ymax && iflag[lw] == 0 ) {
+ if (w[lw][1] >= xmin && w[lw][2] >= ymin && w[lw][1]
+ <= xmax && w[lw][2] <= ymax && iflag[lw] == 0) {
- lwout = lw / ldemo;
- lwout *= ldemo;
- if(lwout == lw && (i == miter || i == iter1)) { /* check this once again */
- stack[nstack][1] = mixx / conv + w[lw][1] / conv;
- stack[nstack][2] = miyy /conv + w[lw][2] / conv;
- stack[nstack][3] = w[lw][3];
+ lwout = lw / ldemo;
+ lwout *= ldemo;
+ if (lwout == lw && (i == miter || i == iter1)) { /* check this once again */
+ stack[nstack][1] = mixx / conv + w[lw][1] / conv;
+ stack[nstack][2] = miyy / conv + w[lw][2] / conv;
+ stack[nstack][3] = w[lw][3];
- nstack++;
-/* iflag[lw] = 1;*/
- }
+ nstack++;
+ /* iflag[lw] = 1; */
+ }
}
-/* walker is leaving the box */
+ /* walker is leaving the box */
-/* if (w[lw][1] <= mixx && w[lw][2] <= bymi && w[lw][1]
- >= bxma && w[lw][2] >= byma && iflag[lw] == 1) {
- iflag[lw] = 0;
- } */
+ /* if (w[lw][1] <= mixx && w[lw][2] <= bymi && w[lw][1]
+ >= bxma && w[lw][2] >= byma && iflag[lw] == 1) {
+ iflag[lw] = 0;
+ } */
-/* every iterout iteration write out the selected ldemo walkers */
-/* and water depth and time */
+ /* every iterout iteration write out the selected ldemo walkers */
+ /* and water depth and time */
- } /* lw loop */
+ } /* lw loop */
if (i == iter1) {
-/* call output for iteration output */
+ /* call output for iteration output */
if (ts == 1) {
- if (erdep != NULL)
- erod(gama); /* divergence of gama field */
+ if (erdep != NULL)
+ erod(gama); /* divergence of gama field */
- conn = (double) nblock / (double) iblock;
- itime = (int) (i * deltap * timec);
- ii=output_data (itime,conn);
- nstack = 0;
- if (ii != 1)
- G_fatal_error (_("Unable to write raster maps"));
+ conn = (double)nblock / (double)iblock;
+ itime = (int)(i * deltap * timec);
+ ii = output_data(itime, conn);
+ nstack = 0;
+ if (ii != 1)
+ G_fatal_error(_("Unable to write raster maps"));
}
}
- if (sfile != NULL) { /* ascii data site file output for gamma - hydrograph*/
+ if (sfile != NULL) { /* ascii data site file output for gamma - hydrograph */
- for (p=0; p < npoints;p++) {
+ for (p = 0; p < npoints; p++) {
- l = (int) ((points[p].east - mixx + stxm) / stepx) - mx -1;
- k = (int) ((points[p].north - miyy + stym) / stepy) - my -1;
+ l = (int)((points[p].east - mixx + stxm) / stepx) - mx -
+ 1;
+ k = (int)((points[p].north - miyy + stym) / stepy) - my -
+ 1;
- if (zz[k][l] != UNDEF) {
+ if (zz[k][l] != UNDEF) {
- if(wdepth == NULL) {
- points[p].z1 = step * gama[k][l] * cchez[k][l]; /* cchez incl. sqrt(sinsl) */
- } else
- points[p].z1 = gama[k][l] * slope[k][l]; /* sediment */
+ if (wdepth == NULL) {
+ points[p].z1 = step * gama[k][l] * cchez[k][l]; /* cchez incl. sqrt(sinsl) */
+ }
+ else
+ points[p].z1 = gama[k][l] * slope[k][l]; /* sediment */
- G_debug(2, " k,l,z1 %d %d %f",k,l,points[p].z1);
+ G_debug(2, " k,l,z1 %d %d %f", k, l, points[p].z1);
- fprintf(fw,"%f %f %f\n",points[p].east/conv,points[p].north/conv,points[p].z1);
- } /*defined area*/
+ fprintf(fw, "%f %f %f\n", points[p].east / conv,
+ points[p].north / conv, points[p].z1);
+ } /*defined area */
- }
+ }
- }
-
+ }
- } /* miter */
-L_800:
-/* if (iwrib != nblock) {
- icount = icoub / iwrib;
+ } /* miter */
- if (icoub == (icount * iwrib)) {
- ++icfl;
- nflw = icfl + 50;
- conn = (double) nblock / (double) iblock;
+ L_800:
+ /* if (iwrib != nblock) {
+ icount = icoub / iwrib;
- }
- }*/
+ if (icoub == (icount * iwrib)) {
+ ++icfl;
+ nflw = icfl + 50;
+ conn = (double) nblock / (double) iblock;
+ }
+ } */
-if (err != NULL)
- {
- for (k = 0; k < my; k++) {
- for (l = 0; l < mx; l++) {
+
+ if (err != NULL) {
+ for (k = 0; k < my; k++) {
+ for (l = 0; l < mx; l++) {
if (zz[k][l] != UNDEF) {
- d1 = gama[k][l] * (double) conn;
- gammas[k][l] += pow(d1, 3./5.);
- } /* DEFined area */
- }
- }
- }
- if (erdep != NULL)
- erod(gama);
- }
-/* ........ end of iblock loop */
+ d1 = gama[k][l] * (double)conn;
+ gammas[k][l] += pow(d1, 3. / 5.);
+ } /* DEFined area */
+ }
+ }
+ }
+ if (erdep != NULL)
+ erod(gama);
+ }
+ /* ........ end of iblock loop */
-}
+}
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c 2008-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/input.c 2008-09-13 08:32:00 UTC (rev 33426)
@@ -12,7 +12,7 @@
/* ************************************************************** */
-/* GRASS input procedures, allocations */
+/* GRASS input procedures, allocations */
/* *************************************************************** */
/*!
@@ -22,729 +22,750 @@
*/
-int input_data (void)
+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. */
- if (sfile != NULL) {
- /* open ascii file for ascii output of gamma */
- fw=fopen("simwe_data.txt","w");
+ /* 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");
- 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();
+ /* 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();
- /* 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||rain_val >= 0.0){
- si = (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 || 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));
+ 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(wdepth!=NULL)
- gama = (double **)G_malloc (sizeof(double *)*(my));
+ if (tauin != NULL)
+ tau = (float **)G_malloc(sizeof(float *) * (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));
+ if (wdepth != NULL)
+ gama = (double **)G_malloc(sizeof(double *) * (my));
- if(rain != NULL||rain_val >= 0.0)
- si[l] = (double*)G_malloc (sizeof(double)*(mx));
+ 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));
- if(infil != NULL||infil_val >= 0.0)
- inf[l] = (double*)G_malloc (sizeof(double)*(mx));
+ if (rain != NULL || rain_val >= 0.0)
+ si[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 (infil != NULL || infil_val >= 0.0)
+ inf[l] = (double *)G_malloc(sizeof(double) * (mx));
- if(detin!=NULL)
- dc[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(tranin!=NULL)
- ct[l] = (float*)G_malloc (sizeof(float)*(mx));
+ if (detin != NULL)
+ dc[l] = (float *)G_malloc(sizeof(float) * (mx));
- if(tauin!=NULL)
- tau[l] = (float*)G_malloc (sizeof(float)*(mx));
+ if (tranin != NULL)
+ ct[l] = (float *)G_malloc(sizeof(float) * (mx));
- if(wdepth!=NULL)
- gama[l] = (double*)G_malloc (sizeof(double)*(mx));
- }
+ if (tauin != NULL)
+ tau[l] = (float *)G_malloc(sizeof(float) * (mx));
- G_debug(3, "Running MAY 10 version, started modifications on 20080211");
+ if (wdepth != NULL)
+ gama[l] = (double *)G_malloc(sizeof(double) * (mx));
+ }
- /* 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);
+ G_debug(3, "Running MAY 10 version, started modifications on 20080211");
- fd1 = G_open_cell_old(elevin,mapset);
+ /* 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);
- /* 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);
+ fd1 = G_open_cell_old(elevin, mapset);
- fd2 = G_open_cell_old(dxin,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(dyin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), dyin);
+ fd2 = G_open_cell_old(dxin, mapset);
- fd3 = G_open_cell_old(dyin,mapset);
- /* END OF REPLACEMENT */
+ if ((mapset = G_find_cell(dyin, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), dyin);
- /* 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);
- }
+ fd3 = G_open_cell_old(dyin, mapset);
+ /* END OF REPLACEMENT */
- /* 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);
- }
+ /* 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(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);
- }
+ /* 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);
+ }
- 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);
+ fd4a = G_open_cell_old(infil, mapset);
+ }
- fd4b = G_open_cell_old(traps,mapset);
- }
+ if (traps != NULL) {
+ if ((mapset = G_find_cell(traps, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), traps);
- if(detin != NULL){
- if((mapset=G_find_cell(detin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), detin);
+ fd4b = G_open_cell_old(traps, mapset);
+ }
- fd9 = G_open_cell_old(detin,mapset);
- }
+ if (detin != NULL) {
+ if ((mapset = G_find_cell(detin, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), detin);
- if(tranin != NULL){
- if((mapset=G_find_cell(tranin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), tranin);
+ fd9 = G_open_cell_old(detin, mapset);
+ }
- fd10 = G_open_cell_old(tranin,mapset);
- }
+ if (tranin != NULL) {
+ if ((mapset = G_find_cell(tranin, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), tranin);
- if(tauin != NULL){
- if((mapset=G_find_cell(tauin,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), tauin);
+ fd10 = G_open_cell_old(tranin, mapset);
+ }
- fd11 = G_open_cell_old(tauin,mapset);
- }
+ if (tauin != NULL) {
+ if ((mapset = G_find_cell(tauin, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), tauin);
- if(wdepth != NULL){
- if((mapset=G_find_cell(wdepth,""))==NULL)
- G_fatal_error(_("Raster map <%s> not found"), wdepth);
+ fd11 = G_open_cell_old(tauin, mapset);
+ }
- fd12 = G_open_cell_old(wdepth,mapset);
- }
+ if (wdepth != NULL) {
+ if ((mapset = G_find_cell(wdepth, "")) == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), wdepth);
- 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);
+ fd12 = G_open_cell_old(wdepth, mapset);
+ }
- if(manin != NULL)
- G_get_f_raster_row(fd5,cell5,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(rain != NULL)
- G_get_d_raster_row(fd4,cell4,row);
+ if (manin != NULL)
+ G_get_f_raster_row(fd5, cell5, 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);
- 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 elevation data exists store in zz[][]*/
- if(!G_is_f_null_value(cell1+j))
- zz[row_rev][j] = (float ) (conv * cell1[j]);
- else
- zz[row_rev][j] = UNDEF;
+ if (wdepth != NULL)
+ G_get_d_raster_row(fd12, cell12, row);
- if(!G_is_d_null_value(cell2+j))
- v1[row_rev][j] = (double ) cell2[j];
- else
- v1[row_rev][j] = UNDEF;
+ 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
+ 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;
- /* undef all area if something's missing */
- if(v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF)
- 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?*/
+ if (!G_is_d_null_value(cell3 + j))
+ v2[row_rev][j] = (double)cell3[j];
+ else
+ v2[row_rev][j] = UNDEF;
- /* 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*/
- /*printf("\n INPUTrain, convert %f %f",si[row_rev][j],unitconv); */
+ /* undef all area if something's missing */
+ if (v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF)
+ zz[row_rev][j] = UNDEF;
- 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*/
- /*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;
- }
- }
+ /* 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 (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 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 */
+ /*printf("\n INPUTrain, convert %f %f",si[row_rev][j],unitconv); */
- 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;
- }
- }
+ else {
+ si[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;
- }
- }
+ /* 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 */
+ /*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 (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;
- }
- }
+ 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 (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 (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 (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);
+ }
+ G_close_cell(fd1);
+ G_close_cell(fd2);
+ G_close_cell(fd3);
- if(rain != NULL)
- G_close_cell(fd4);
+ if (rain != NULL)
+ G_close_cell(fd4);
- if(infil != NULL)
- G_close_cell(fd4a);
+ if (infil != NULL)
+ G_close_cell(fd4a);
- if(traps != NULL)
- G_close_cell(fd4b);
- /* Maybe a conditional to manin!=NULL here ! */
- G_close_cell(fd5);
+ if (traps != NULL)
+ G_close_cell(fd4b);
+ /* Maybe a conditional to manin!=NULL here ! */
+ G_close_cell(fd5);
+
/****************/
- if(detin != NULL)
- G_close_cell(fd9);
+ if (detin != NULL)
+ G_close_cell(fd9);
- if(tranin != NULL)
- G_close_cell(fd10);
+ if (tranin != NULL)
+ G_close_cell(fd10);
- if(tauin != NULL)
- G_close_cell(fd11);
+ if (tauin != NULL)
+ G_close_cell(fd11);
- if(wdepth != NULL)
- G_close_cell(fd12);
+ if (wdepth != NULL)
+ G_close_cell(fd12);
- return 1;
+ return 1;
}
/* data preparations, sigma, shear, etc. */
-int grad_check (void)
+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;
-
- sisum = 0.;
- infsum = 0.;
- cmul2 = rhow * gacc;
-
- /* mandatory alloc. - should be moved to main.c*/
- slope = (double **)G_malloc (sizeof(double *)*(my));
+ 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;
- for(l=0;l<my;l++)
- slope[l] = (double*)G_malloc (sizeof(double)*(mx));
+ sisum = 0.;
+ infsum = 0.;
+ cmul2 = rhow * gacc;
- for (j = 0; j < my; j++)
- {
- for (i = 0; i < mx; i++)
- slope[j][i] = 0.;
- }
+ /* 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 (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.);
- /* 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 */
+ 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 (inf != NULL && smax < infmax)
- G_warning (_("Infiltration exceeds the rainfall rate everywhere! No overland flow."));
+ }
+ if (inf != NULL && smax < infmax)
+ G_warning(_("Infiltration exceeds the rainfall rate everywhere! No overland flow."));
- cc = (double) mx * my;
- si0 = sisum / cc;
- vmean = vsum / cc;
- chmean = chsum / cc;
+ cc = (double)mx *my;
- if (inf != NULL) infmean = infsum / cc;
+ si0 = sisum / cc;
+ vmean = vsum / cc;
+ chmean = chsum / cc;
- if (wdepth != NULL) deltaw = 0.8/(sigmax * vmax); /*time step for sediment*/
- deltap = 0.25 * sqrt(stepx * stepy)/vmean; /*time step for water*/
+ if (inf != NULL)
+ infmean = infsum / cc;
- if(deltaw > deltap)
- timec = 4.;
- else
- timec = 1.25;
+ if (wdepth != NULL)
+ deltaw = 0.8 / (sigmax * vmax); /*time step for sediment */
+ deltap = 0.25 * sqrt(stepx * stepy) / vmean; /*time step for water */
- 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 Source Rate (rainf. excess or sediment) \t= %f m/s or kg/m2s \n"),si0);
- G_message (_("Mean flow velocity \t= %f m/s\n"),vmean);
- G_message (_("Mean Mannings \t= %f\n"),1.0/chmean);
+ if (deltaw > deltap)
+ timec = 4.;
+ else
+ timec = 1.25;
- deltap = amin1(deltap,deltaw);
+ 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 */
- 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 */
+ fprintf(stderr, "\n");
+ G_message(_("Min elevation \t= %.2f m\nMax elevation \t= %.2f m\n"), zmin,
+ zmax);
+ G_message(_("Mean Source Rate (rainf. excess or sediment) \t= %f m/s or kg/m2s \n"),
+ si0);
+ G_message(_("Mean flow velocity \t= %f m/s\n"), vmean);
+ G_message(_("Mean Mannings \t= %f\n"), 1.0 / chmean);
+ deltap = amin1(deltap, deltaw);
- /*! 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 */
+ 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 */
+
+
+ /*! 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 */
}
+ }
- /*! 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"));
- }
+ /*! 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"));
+ }
- /*! 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) {
- /* 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 */
+ /*! 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) {
+ /* 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 */
}
- return 1;
+ }
+ return 1;
}
-double amax1 (double arg1, double arg2)
+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 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 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 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-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/output.c 2008-09-13 08:32:00 UTC (rev 33426)
@@ -12,707 +12,740 @@
int output_data(int tt, double ft)
-
{
- FCELL *cell6, *cell7, *cell8;
- FCELL *cell14, *cell15, *cell16;
- int fd6, fd7, fd8;
- int fd14, fd15,fd16;
- int i,iarc,j;
- float gsmax=0,dismax=0.,gmax=0.,ermax=-1.e+12,ermin=1.e+12;
- struct Colors colors;
- struct History hist, hist1; /* hist2, hist3, hist4, hist5*/
- char *depth0=NULL,*disch0=NULL,*err0=NULL;
- char *conc0=NULL,*flux0=NULL;
- char *erdep0=NULL,*outwalk0=NULL;
- char *mapst = NULL;
- char *type;
- char buf[256];
- int ndigit;
- FCELL dat1, dat2;
- float a1,a2;
- Site_head walkershead;
- Site *sd;
+ FCELL *cell6, *cell7, *cell8;
+ FCELL *cell14, *cell15, *cell16;
+ int fd6, fd7, fd8;
+ int fd14, fd15, fd16;
+ int i, iarc, j;
+ float gsmax = 0, dismax = 0., gmax = 0., ermax = -1.e+12, ermin = 1.e+12;
+ struct Colors colors;
+ struct History hist, hist1; /* hist2, hist3, hist4, hist5 */
+ char *depth0 = NULL, *disch0 = NULL, *err0 = NULL;
+ char *conc0 = NULL, *flux0 = NULL;
+ char *erdep0 = NULL, *outwalk0 = NULL;
+ char *mapst = NULL;
+ char *type;
+ char buf[256];
+ int ndigit;
+ FCELL dat1, dat2;
+ float a1, a2;
+ Site_head walkershead;
+ Site *sd;
- ndigit=2;
- if(timesec >= 10) ndigit=3;
- if(timesec >= 100) ndigit=4;
- if(timesec >= 1000) ndigit=5;
- if(timesec >= 10000) ndigit=6;
+ ndigit = 2;
+ if (timesec >= 10)
+ ndigit = 3;
+ if (timesec >= 100)
+ ndigit = 4;
+ if (timesec >= 1000)
+ ndigit = 5;
+ if (timesec >= 10000)
+ ndigit = 6;
- if (outwalk != NULL)
- {
- if (ts == 1) {
- sprintf(buf,"%s%.*d",outwalk,ndigit,tt);
- outwalk0 = G_store(buf);
- fdoutwalk = G_fopen_sites_new (outwalk0);
- }
- else
- fdoutwalk = G_fopen_sites_new (outwalk);
+ if (outwalk != NULL) {
+ if (ts == 1) {
+ sprintf(buf, "%s%.*d", outwalk, ndigit, tt);
+ outwalk0 = G_store(buf);
+ fdoutwalk = G_fopen_sites_new(outwalk0);
+ }
+ else
+ fdoutwalk = G_fopen_sites_new(outwalk);
- if (fdoutwalk == NULL)
- G_fatal_error ("Cannot open %s", outwalk);
- else
- {
- if(NULL == (sd = G_site_new_struct(-1,2,0,1)))
- G_fatal_error("memory allocation failed for site");
+ if (fdoutwalk == NULL)
+ G_fatal_error("Cannot open %s", outwalk);
+ else {
+ if (NULL == (sd = G_site_new_struct(-1, 2, 0, 1)))
+ G_fatal_error("memory allocation failed for site");
- if (ts == 1)
- walkershead.name = outwalk0;
- else
- walkershead.name = outwalk;
+ if (ts == 1)
+ walkershead.name = outwalk0;
+ else
+ walkershead.name = outwalk;
- walkershead.desc = G_strdup ("output walkers");
- walkershead.desc = (char *) G_malloc (128 * sizeof (char));
- sprintf (walkershead.desc, "output walkers of %s [raster]",
- depth);
- walkershead.time = NULL;
- walkershead.stime = NULL;
- walkershead.labels = NULL;
- walkershead.form = NULL;
+ walkershead.desc = G_strdup("output walkers");
+ walkershead.desc = (char *)G_malloc(128 * sizeof(char));
+ sprintf(walkershead.desc, "output walkers of %s [raster]", depth);
+ walkershead.time = NULL;
+ walkershead.stime = NULL;
+ walkershead.labels = NULL;
+ walkershead.form = NULL;
- G_site_put_head (fdoutwalk, &walkershead);
+ G_site_put_head(fdoutwalk, &walkershead);
- for (i = 0; i < nstack; i++) {
- sd->east = (float)stack[i][1];
- sd->north = (float)stack[i][2];
- sd->fcat = (float)stack[i][3];
- G_site_put(fdoutwalk, sd);
+ for (i = 0; i < nstack; i++) {
+ sd->east = (float)stack[i][1];
+ sd->north = (float)stack[i][2];
+ sd->fcat = (float)stack[i][3];
+ G_site_put(fdoutwalk, sd);
+ }
+
}
-
}
- }
- if (depth != NULL)
- {
- cell6 = G_allocate_f_raster_buf();
- if (ts == 1) {
- sprintf(buf,"%s.%.*d",depth,ndigit,tt);
- depth0 = G_store(buf);
- fd6 = G_open_fp_cell_new (depth0);
+ if (depth != NULL) {
+ cell6 = G_allocate_f_raster_buf();
+ if (ts == 1) {
+ sprintf(buf, "%s.%.*d", depth, ndigit, tt);
+ depth0 = G_store(buf);
+ fd6 = G_open_fp_cell_new(depth0);
}
else
- fd6 = G_open_fp_cell_new (depth);
- if (fd6 < 0)
- G_fatal_error ("unable to create raster map %s", depth);
- }
+ fd6 = G_open_fp_cell_new(depth);
+ if (fd6 < 0)
+ G_fatal_error("unable to create raster map %s", depth);
+ }
- if (disch != NULL)
- {
- cell7 = G_allocate_f_raster_buf();
- if (ts == 1) {
- sprintf(buf,"%s.%.*d",disch,ndigit,tt);
- disch0 = G_store(buf);
- fd7 = G_open_fp_cell_new (disch0);
+ if (disch != NULL) {
+ cell7 = G_allocate_f_raster_buf();
+ if (ts == 1) {
+ sprintf(buf, "%s.%.*d", disch, ndigit, tt);
+ disch0 = G_store(buf);
+ fd7 = G_open_fp_cell_new(disch0);
}
else
- fd7 = G_open_fp_cell_new (disch);
- if (fd7 < 0)
- G_fatal_error("unable to create raster map %s", disch);
- }
+ fd7 = G_open_fp_cell_new(disch);
+ if (fd7 < 0)
+ G_fatal_error("unable to create raster map %s", disch);
+ }
- if (err != NULL)
- {
- cell8 = G_allocate_f_raster_buf();
- if (ts == 1) {
- sprintf(buf,"%s.%.*d",err,ndigit,tt);
- err0 = G_store(buf);
- fd8 = G_open_fp_cell_new (err0);
+ if (err != NULL) {
+ cell8 = G_allocate_f_raster_buf();
+ if (ts == 1) {
+ sprintf(buf, "%s.%.*d", err, ndigit, tt);
+ err0 = G_store(buf);
+ fd8 = G_open_fp_cell_new(err0);
}
else
- fd8 = G_open_fp_cell_new (err);
+ fd8 = G_open_fp_cell_new(err);
- if (fd8 < 0)
- G_fatal_error ("unable to create raster map %s", err);
- }
+ if (fd8 < 0)
+ G_fatal_error("unable to create raster map %s", err);
+ }
- if (conc != NULL)
- {
- cell14 = G_allocate_f_raster_buf();
- if (ts == 1) {
- sprintf(buf,"%s.%.*d",conc,ndigit,tt);
- conc0 = G_store(buf);
- fd14 = G_open_fp_cell_new (conc0);
- }
- else
- fd14 = G_open_fp_cell_new (conc);
+ if (conc != NULL) {
+ cell14 = G_allocate_f_raster_buf();
+ if (ts == 1) {
+ sprintf(buf, "%s.%.*d", conc, ndigit, tt);
+ conc0 = G_store(buf);
+ fd14 = G_open_fp_cell_new(conc0);
+ }
+ else
+ fd14 = G_open_fp_cell_new(conc);
- if (fd14 < 0)
- G_fatal_error ("unable to create raster map %s", conc);
- }
+ if (fd14 < 0)
+ G_fatal_error("unable to create raster map %s", conc);
+ }
- if (flux != NULL)
- {
- cell15 = G_allocate_f_raster_buf();
- if (ts == 1) {
- sprintf(buf,"%s.%.*d",flux,ndigit,tt);
- flux0 = G_store(buf);
- fd15 = G_open_fp_cell_new (flux0);
- }
- else
- fd15 = G_open_fp_cell_new (flux);
+ if (flux != NULL) {
+ cell15 = G_allocate_f_raster_buf();
+ if (ts == 1) {
+ sprintf(buf, "%s.%.*d", flux, ndigit, tt);
+ flux0 = G_store(buf);
+ fd15 = G_open_fp_cell_new(flux0);
+ }
+ else
+ fd15 = G_open_fp_cell_new(flux);
- if (fd15 < 0)
- G_fatal_error ("unable to create raster map %s", flux);
- }
+ if (fd15 < 0)
+ G_fatal_error("unable to create raster map %s", flux);
+ }
- if (erdep != NULL)
- {
- cell16 = G_allocate_f_raster_buf();
- if (ts == 1) {
- sprintf(buf,"%s.%.*d",erdep,ndigit,tt);
- erdep0 = G_store(buf);
- fd16 = G_open_fp_cell_new (erdep0);
- }
- else
- fd16 = G_open_fp_cell_new (erdep);
+ if (erdep != NULL) {
+ cell16 = G_allocate_f_raster_buf();
+ if (ts == 1) {
+ sprintf(buf, "%s.%.*d", erdep, ndigit, tt);
+ erdep0 = G_store(buf);
+ fd16 = G_open_fp_cell_new(erdep0);
+ }
+ else
+ fd16 = G_open_fp_cell_new(erdep);
- if (fd16 < 0)
- G_fatal_error ("unable to create raster map %s", erdep);
- }
+ if (fd16 < 0)
+ G_fatal_error("unable to create raster map %s", erdep);
+ }
- if(G_set_window (&cellhd) < 0)
- exit(3);
+ if (G_set_window(&cellhd) < 0)
+ exit(3);
- if (my != G_window_rows())
- G_fatal_error ("OOPS: rows changed from %d to %d\n", mx, G_window_rows());
- if (mx != G_window_cols())
- G_fatal_error ("OOPS: cols changed from %d to %d\n", my, G_window_cols());
+ if (my != G_window_rows())
+ G_fatal_error("OOPS: rows changed from %d to %d\n", mx,
+ G_window_rows());
+ if (mx != G_window_cols())
+ G_fatal_error("OOPS: cols changed from %d to %d\n", my,
+ G_window_cols());
- for(iarc=0;iarc<my;iarc++)
- {
- i=my-iarc-1;
- if(depth!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || gama[i][j] == UNDEF)
- G_set_f_null_value (cell6+j,1);
- else {
- a1 = pow(gama[i][j],3./5.);
- cell6[j]=(FCELL) a1;/* add conv? */
- gmax = amax1(gmax, a1);
- }
- }
- G_put_f_raster_row (fd6, cell6);
- }
+ for (iarc = 0; iarc < my; iarc++) {
+ i = my - iarc - 1;
+ if (depth != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF)
+ G_set_f_null_value(cell6 + j, 1);
+ else {
+ a1 = pow(gama[i][j], 3. / 5.);
+ cell6[j] = (FCELL) a1; /* add conv? */
+ gmax = amax1(gmax, a1);
+ }
+ }
+ G_put_f_raster_row(fd6, cell6);
+ }
- if(disch!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || gama[i][j] == UNDEF || cchez[i][j] == UNDEF)
- G_set_f_null_value (cell7+j,1);
- else {
- a2 = step * gama[i][j] * cchez[i][j]; /* cchez incl. sqrt(sinsl) */
- cell7[j]=(FCELL) a2;/* add conv? */
- dismax = amax1(dismax, a2);
- }
- }
- G_put_f_raster_row (fd7, cell7);
- }
+ if (disch != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF ||
+ cchez[i][j] == UNDEF)
+ G_set_f_null_value(cell7 + j, 1);
+ else {
+ a2 = step * gama[i][j] * cchez[i][j]; /* cchez incl. sqrt(sinsl) */
+ cell7[j] = (FCELL) a2; /* add conv? */
+ dismax = amax1(dismax, a2);
+ }
+ }
+ G_put_f_raster_row(fd7, cell7);
+ }
- if(err!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || gammas[i][j] == UNDEF)
- G_set_f_null_value (cell8+j,1);
+ if (err != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gammas[i][j] == UNDEF)
+ G_set_f_null_value(cell8 + j, 1);
else {
- cell8[j]=(FCELL)gammas[i][j];
- gsmax = amax1(gsmax, gammas[i][j]);/* add conv? */
+ cell8[j] = (FCELL) gammas[i][j];
+ gsmax = amax1(gsmax, gammas[i][j]); /* add conv? */
}
- }
- G_put_f_raster_row (fd8, cell8);
- }
+ }
+ G_put_f_raster_row(fd8, cell8);
+ }
- if(conc!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || gama[i][j] == UNDEF)
- G_set_f_null_value (cell14+j,1);
+ if (conc != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF)
+ G_set_f_null_value(cell14 + j, 1);
else {
- cell14[j]=(FCELL)gama[i][j];
- /* gsmax = amax1(gsmax, gama[i][j]); */
+ cell14[j] = (FCELL) gama[i][j];
+ /* gsmax = amax1(gsmax, gama[i][j]); */
}
- }
- G_put_f_raster_row (fd14, cell14);
- }
+ }
+ G_put_f_raster_row(fd14, cell14);
+ }
- if(flux!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || gama[i][j] == UNDEF || slope[i][j] == UNDEF)
- G_set_f_null_value (cell15+j,1);
+ if (flux != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || gama[i][j] == UNDEF ||
+ slope[i][j] == UNDEF)
+ G_set_f_null_value(cell15 + j, 1);
else {
- a2 = gama[i][j] * slope[i][j];
- cell15[j]=(FCELL) a2;
- dismax = amax1(dismax, a2);
- }
+ a2 = gama[i][j] * slope[i][j];
+ cell15[j] = (FCELL) a2;
+ dismax = amax1(dismax, a2);
}
- G_put_f_raster_row (fd15, cell15);
- }
+ }
+ G_put_f_raster_row(fd15, cell15);
+ }
- if(erdep!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || er[i][j] == UNDEF)
- G_set_f_null_value (cell16+j,1);
+ if (erdep != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || er[i][j] == UNDEF)
+ G_set_f_null_value(cell16 + j, 1);
else {
- cell16[j]=(FCELL)er[i][j];
- ermax = amax1(ermax, er[i][j]);
- ermin = amin1(ermin, er[i][j]);
- }
+ cell16[j] = (FCELL) er[i][j];
+ ermax = amax1(ermax, er[i][j]);
+ ermin = amin1(ermin, er[i][j]);
}
- G_put_f_raster_row (fd16, cell16);
- }
-
+ }
+ G_put_f_raster_row(fd16, cell16);
}
- if(depth!=NULL)
- G_close_cell (fd6);
- if(disch!=NULL)
- G_close_cell (fd7);
- if(err!=NULL)
- G_close_cell (fd8);
- if(conc!=NULL)
- G_close_cell (fd14);
- if(flux!=NULL)
- G_close_cell (fd15);
- if(erdep!=NULL)
- G_close_cell (fd16);
+ }
- if (depth!=NULL) {
+ if (depth != NULL)
+ G_close_cell(fd6);
+ if (disch != NULL)
+ G_close_cell(fd7);
+ if (err != NULL)
+ G_close_cell(fd8);
+ if (conc != NULL)
+ G_close_cell(fd14);
+ if (flux != NULL)
+ G_close_cell(fd15);
+ if (erdep != NULL)
+ G_close_cell(fd16);
- G_init_colors(&colors);
+ if (depth != NULL) {
- dat1 = (FCELL) 0.;
- dat2 = (FCELL) 0.001;
- G_add_f_raster_color_rule(&dat1,255,255,255,&dat2,255,255,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.05;
- G_add_f_raster_color_rule(&dat1,255,255,0,&dat2,0,255,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.1;
- G_add_f_raster_color_rule(&dat1,0,255,255,&dat2,0,127,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.5;
- G_add_f_raster_color_rule(&dat1,0,127,255,&dat2,0,0,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) gmax;
- G_add_f_raster_color_rule(&dat1,0,0,255,&dat2,0,0,0,&colors);
+ G_init_colors(&colors);
+ dat1 = (FCELL) 0.;
+ dat2 = (FCELL) 0.001;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.05;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 0, &dat2, 0, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.1;
+ G_add_f_raster_color_rule(&dat1, 0, 255, 255, &dat2, 0, 127, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.5;
+ G_add_f_raster_color_rule(&dat1, 0, 127, 255, &dat2, 0, 0, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) gmax;
+ G_add_f_raster_color_rule(&dat1, 0, 0, 255, &dat2, 0, 0, 0, &colors);
- if (ts == 1) {
+
+ if (ts == 1) {
if ((mapst = G_find_file("fcell", depth0, "")) == NULL)
- G_fatal_error("cannot find file %s", depth0);
+ G_fatal_error("cannot find file %s", depth0);
G_write_colors(depth0, mapst, &colors);
- G_quantize_fp_map_range(depth0,mapst,0.,(FCELL)gmax,0,(CELL)gmax);
+ G_quantize_fp_map_range(depth0, mapst, 0., (FCELL) gmax, 0,
+ (CELL) gmax);
G_free_colors(&colors);
- }
+ }
else {
if ((mapst = G_find_file("fcell", depth, "")) == NULL)
- G_fatal_error("cannot find file %s", depth);
+ G_fatal_error("cannot find file %s", depth);
G_write_colors(depth, mapst, &colors);
- G_quantize_fp_map_range(depth,mapst,0.,(FCELL)gmax,0,(CELL)gmax);
+ G_quantize_fp_map_range(depth, mapst, 0., (FCELL) gmax, 0,
+ (CELL) gmax);
G_free_colors(&colors);
- }
+ }
- }
+ }
- if (disch!=NULL) {
+ if (disch != NULL) {
- G_init_colors(&colors);
+ G_init_colors(&colors);
- dat1 = (FCELL) 0.;
- dat2 = (FCELL) 0.0005;
- G_add_f_raster_color_rule(&dat1,255,255,255,&dat2,255,255,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.005;
- G_add_f_raster_color_rule(&dat1,255,255,0,&dat2,0,255,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.05;
- G_add_f_raster_color_rule(&dat1,0,255,255,&dat2,0,127,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.1;
- G_add_f_raster_color_rule(&dat1,0,127,255,&dat2,0,0,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) dismax;
- G_add_f_raster_color_rule(&dat1,0,0,255,&dat2,0,0,0,&colors);
+ dat1 = (FCELL) 0.;
+ dat2 = (FCELL) 0.0005;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.005;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 0, &dat2, 0, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.05;
+ G_add_f_raster_color_rule(&dat1, 0, 255, 255, &dat2, 0, 127, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.1;
+ G_add_f_raster_color_rule(&dat1, 0, 127, 255, &dat2, 0, 0, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) dismax;
+ G_add_f_raster_color_rule(&dat1, 0, 0, 255, &dat2, 0, 0, 0, &colors);
- if (ts == 1) {
+ if (ts == 1) {
if ((mapst = G_find_file("cell", disch0, "")) == NULL)
- G_fatal_error ("cannot find file %s", disch0);
+ G_fatal_error("cannot find file %s", disch0);
G_write_colors(disch0, mapst, &colors);
- G_quantize_fp_map_range(disch0,mapst,0.,(FCELL)dismax,0,(CELL)dismax);
- G_free_colors(&colors);
- }
+ G_quantize_fp_map_range(disch0, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
+ G_free_colors(&colors);
+ }
else {
if ((mapst = G_find_file("cell", disch, "")) == NULL)
- G_fatal_error ("cannot find file %s", disch);
+ G_fatal_error("cannot find file %s", disch);
G_write_colors(disch, mapst, &colors);
- G_quantize_fp_map_range(disch, mapst,0.,(FCELL)dismax,0,(CELL)dismax);
+ G_quantize_fp_map_range(disch, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
G_free_colors(&colors);
- }
- }
+ }
+ }
- if (flux!=NULL) {
+ if (flux != NULL) {
- G_init_colors(&colors);
+ G_init_colors(&colors);
- dat1 = (FCELL) 0.;
- dat2 = (FCELL) 0.001;
- G_add_f_raster_color_rule(&dat1,255,255,255,&dat2,255,255,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.1;
- G_add_f_raster_color_rule(&dat1,255,255,0,&dat2,255,127,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 1.;
- G_add_f_raster_color_rule(&dat1,255,127,0,&dat2,191,127,63,&colors);
- dat1 = dat2;
- dat2 = (FCELL) dismax;
- G_add_f_raster_color_rule(&dat1,191,127,63,&dat2,0,0,0,&colors);
+ dat1 = (FCELL) 0.;
+ dat2 = (FCELL) 0.001;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.1;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 1.;
+ G_add_f_raster_color_rule(&dat1, 255, 127, 0, &dat2, 191, 127, 63,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) dismax;
+ G_add_f_raster_color_rule(&dat1, 191, 127, 63, &dat2, 0, 0, 0,
+ &colors);
- if (ts == 1) {
- if ((mapst = G_find_file("cell", flux0, "")) == NULL)
- G_fatal_error ("cannot find file %s", flux0);
+ if (ts == 1) {
+ if ((mapst = G_find_file("cell", flux0, "")) == NULL)
+ G_fatal_error("cannot find file %s", flux0);
G_write_colors(flux0, mapst, &colors);
- G_quantize_fp_map_range(flux0,mapst,0.,(FCELL)dismax,0,(CELL)dismax);
+ G_quantize_fp_map_range(flux0, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
G_free_colors(&colors);
+ }
+ else {
+
+ if ((mapst = G_find_file("cell", flux, "")) == NULL)
+ G_fatal_error("cannot find file %s", flux);
+ G_write_colors(flux, mapst, &colors);
+ G_quantize_fp_map_range(flux, mapst, 0., (FCELL) dismax, 0,
+ (CELL) dismax);
+ G_free_colors(&colors);
+ }
}
- else {
-
- if ((mapst = G_find_file("cell", flux, "")) == NULL)
- G_fatal_error ("cannot find file %s", flux);
- G_write_colors(flux, mapst, &colors);
- G_quantize_fp_map_range(flux, mapst,0.,(FCELL)dismax,0,(CELL)dismax);
- G_free_colors(&colors);
- }
- }
- if (erdep!=NULL) {
+ if (erdep != NULL) {
- G_init_colors(&colors);
+ G_init_colors(&colors);
- dat1 = (FCELL) ermax;
- dat2 = (FCELL) 0.1;
- G_add_f_raster_color_rule(&dat1,0,0,0,&dat2,0,0,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.01;
- G_add_f_raster_color_rule(&dat1,0,0,255,&dat2,0,191,191,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.0001;
- G_add_f_raster_color_rule(&dat1,0,191,191,&dat2,170,255,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.;
- G_add_f_raster_color_rule(&dat1,170,255,255,&dat2,255,255,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) -0.0001;
- G_add_f_raster_color_rule(&dat1,255,255,255,&dat2,255,255,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) -0.01;
- G_add_f_raster_color_rule(&dat1,255,255,0,&dat2,255,127,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) -0.1;
- G_add_f_raster_color_rule(&dat1,255,127,0,&dat2,255,0,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) ermin;
- G_add_f_raster_color_rule(&dat1,255,0,0,&dat2,255,0,255,&colors);
+ dat1 = (FCELL) ermax;
+ dat2 = (FCELL) 0.1;
+ G_add_f_raster_color_rule(&dat1, 0, 0, 0, &dat2, 0, 0, 255, &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.01;
+ G_add_f_raster_color_rule(&dat1, 0, 0, 255, &dat2, 0, 191, 191,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.0001;
+ G_add_f_raster_color_rule(&dat1, 0, 191, 191, &dat2, 170, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.;
+ G_add_f_raster_color_rule(&dat1, 170, 255, 255, &dat2, 255, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.0001;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.01;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.1;
+ G_add_f_raster_color_rule(&dat1, 255, 127, 0, &dat2, 255, 0, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) ermin;
+ G_add_f_raster_color_rule(&dat1, 255, 0, 0, &dat2, 255, 0, 255,
+ &colors);
- if (ts == 1) {
- if ((mapst = G_find_file("cell", erdep0, "")) == NULL)
- G_fatal_error ("cannot find file %s", erdep0);
- G_write_colors(erdep0, mapst, &colors);
- G_quantize_fp_map_range(erdep0,mapst,(FCELL)ermin,(FCELL)ermax,(CELL)ermin,(CELL)ermax);
- G_free_colors(&colors);
+ if (ts == 1) {
+ if ((mapst = G_find_file("cell", erdep0, "")) == NULL)
+ G_fatal_error("cannot find file %s", erdep0);
+ G_write_colors(erdep0, mapst, &colors);
+ G_quantize_fp_map_range(erdep0, mapst, (FCELL) ermin,
+ (FCELL) ermax, (CELL) ermin,
+ (CELL) ermax);
+ G_free_colors(&colors);
type = "raster";
- G_short_history (erdep0, type, &hist1);
- sprintf (hist1.edhist[0], "The sediment flux file is %s", flux0);
- hist1.edlinecnt = 1;
- G_write_history (erdep0, &hist1);
- }
- else {
+ G_short_history(erdep0, type, &hist1);
+ sprintf(hist1.edhist[0], "The sediment flux file is %s", flux0);
+ hist1.edlinecnt = 1;
+ G_write_history(erdep0, &hist1);
+ }
+ else {
- if ((mapst = G_find_file("cell", erdep, "")) == NULL)
- G_fatal_error ("cannot find file %s", erdep);
- G_write_colors(erdep, mapst, &colors);
- G_quantize_fp_map_range(erdep, mapst,(FCELL)ermin,(FCELL)ermax,(CELL)ermin,(CELL)ermax);
- G_free_colors(&colors);
+ if ((mapst = G_find_file("cell", erdep, "")) == NULL)
+ G_fatal_error("cannot find file %s", erdep);
+ G_write_colors(erdep, mapst, &colors);
+ G_quantize_fp_map_range(erdep, mapst, (FCELL) ermin,
+ (FCELL) ermax, (CELL) ermin,
+ (CELL) ermax);
+ G_free_colors(&colors);
- type = "raster";
- G_short_history (erdep, type, &hist1);
- sprintf (hist1.edhist[0], "The sediment flux file is %s", flux);
- hist1.edlinecnt = 1;
- G_write_history (erdep, &hist1);
+ type = "raster";
+ G_short_history(erdep, type, &hist1);
+ sprintf(hist1.edhist[0], "The sediment flux file is %s", flux);
+ hist1.edlinecnt = 1;
+ G_write_history(erdep, &hist1);
+ }
}
- }
-/* history section */
+ /* history section */
- if (depth != NULL)
- {
- type = "raster";
+ if (depth != NULL) {
+ type = "raster";
if (ts == 0) {
- mapst = G_find_file ("cell", depth, "");
- if (mapst == NULL)
- {
- G_warning("File [%s] not found", depth);
- return -1;
- }
- G_short_history (depth, type, &hist);
- } else
- G_short_history (depth0, type, &hist);
+ mapst = G_find_file("cell", depth, "");
+ if (mapst == NULL) {
+ G_warning("File [%s] not found", depth);
+ return -1;
+ }
+ G_short_history(depth, type, &hist);
+ }
+ else
+ G_short_history(depth0, type, &hist);
-/* fprintf (stdout, "\n history initiated\n");
- fflush(stdout);*/
+ /* 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[1],"duration (sec.)=%d, time-serie iteration=%d",timesec,tt);
+ sprintf(hist.edhist[0],
+ "init.walk=%d, maxwalk=%d, remaining 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",
- lwwfin, deltap, vmean);
+ sprintf(hist.edhist[2], "written walkers=%d, deltap=%f, mean vel.=%f",
+ lwwfin, deltap, vmean);
- sprintf(hist.edhist[3], "mean source (si)=%e, mean infil=%e", si0, infmean);
+ sprintf(hist.edhist[3], "mean source (si)=%e, mean infil=%e", si0,
+ infmean);
- 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;
+ 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);
+ G_command_history(&hist);
- if (ts == 1) G_write_history (depth0, &hist);
+ if (ts == 1)
+ G_write_history(depth0, &hist);
else
- G_write_history (depth, &hist);
- }
+ G_write_history(depth, &hist);
+ }
- if (disch != NULL)
- {
- type = "raster";
+ if (disch != NULL) {
+ type = "raster";
if (ts == 0) {
- mapst = G_find_file ("cell", disch, "");
- if (mapst == NULL)
- G_fatal_error ("file [%s] not found\n", disch);
- G_short_history (disch, type, &hist);
- } else
- G_short_history (disch0, type, &hist);
+ mapst = G_find_file("cell", disch, "");
+ if (mapst == NULL)
+ G_fatal_error("file [%s] not found\n", disch);
+ G_short_history(disch, type, &hist);
+ }
+ else
+ G_short_history(disch0, type, &hist);
-/* fprintf (stdout, "\n history initiated\n");
- fflush(stdout);*/
+ /* fprintf (stdout, "\n history initiated\n");
+ fflush(stdout); */
- 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[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",
- lwwfin, deltap, vmean);
+ sprintf(hist.edhist[2], "written walkers=%d, deltap=%f, mean vel.=%f",
+ lwwfin, deltap, vmean);
- sprintf(hist.edhist[3], "mean source (si)=%e, mean infil=%e", si0, infmean);
+ sprintf(hist.edhist[3], "mean source (si)=%e, mean infil=%e", si0,
+ infmean);
- 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;
+ 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);
+ G_command_history(&hist);
- if (ts == 1) G_write_history (disch0, &hist);
+ if (ts == 1)
+ G_write_history(disch0, &hist);
else
- G_write_history (disch, &hist);
- }
+ G_write_history(disch, &hist);
+ }
- if (flux != NULL)
- {
- type = "raster";
- if (ts == 0) {
- mapst = G_find_file ("cell", flux, "");
- if (mapst == NULL)
- G_fatal_error ("file [%s] not found\n", flux);
- G_short_history (flux, type, &hist);
- } else
- G_short_history (flux0, type, &hist);
+ if (flux != NULL) {
+ type = "raster";
+ if (ts == 0) {
+ mapst = G_find_file("cell", flux, "");
+ if (mapst == NULL)
+ G_fatal_error("file [%s] not found\n", flux);
+ G_short_history(flux, type, &hist);
+ }
+ else
+ G_short_history(flux0, type, &hist);
-/* fprintf (stdout, "\n history initiated\n");
- fflush(stdout);*/
+ /* 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[1],"duration (sec.)=%d, time-serie iteration=%d",timesec,tt);
+ sprintf(hist.edhist[0],
+ "init.walk=%d, maxwalk=%d, remaining 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",
- lwwfin, deltap, vmean);
+ sprintf(hist.edhist[2], "written walkers=%d, deltap=%f, mean vel.=%f",
+ lwwfin, deltap, vmean);
- sprintf(hist.edhist[3], "mean source (si)=%f", si0);
+ sprintf(hist.edhist[3], "mean source (si)=%f", si0);
- sprintf (hist.datsrc_1, "input files: %s %s %s", wdepth,dxin,dyin);
- sprintf (hist.datsrc_2, "input files: %s %s %s %s", manin,detin,tranin,tauin);
+ sprintf(hist.datsrc_1, "input files: %s %s %s", wdepth, dxin, dyin);
+ sprintf(hist.datsrc_2, "input files: %s %s %s %s", manin, detin,
+ tranin, tauin);
- hist.edlinecnt = 4;
+ hist.edlinecnt = 4;
- G_command_history(&hist);
+ G_command_history(&hist);
- if (ts == 1) G_write_history (flux0, &hist);
+ if (ts == 1)
+ G_write_history(flux0, &hist);
else
- G_write_history (flux, &hist);
- }
+ G_write_history(flux, &hist);
+ }
- return 1;
+ return 1;
}
int output_et()
-
{
- FCELL *cell13,*cell17;
- int fd13,fd17;
- int i,iarc,j;
- float etmax=-1.e+12,etmin=1.e+12;
- float trc;
- struct Colors colors;
- char *mapst = NULL;
-/* char buf[256]; */
- FCELL dat1, dat2;
-/* float a1,a2; */
+ FCELL *cell13, *cell17;
+ int fd13, fd17;
+ int i, iarc, j;
+ float etmax = -1.e+12, etmin = 1.e+12;
+ float trc;
+ struct Colors colors;
+ char *mapst = NULL;
+ /* char buf[256]; */
+ FCELL dat1, dat2;
+ /* float a1,a2; */
- if (et != NULL)
- {
- cell17 = G_allocate_f_raster_buf();
-/* if (ts == 1) {
- sprintf(buf,"%s.%.*d",et,ndigit,tt);
- et0 = G_store(buf);
- fd17 = G_open_fp_cell_new (et0);
- }
- else*/
- fd17 = G_open_fp_cell_new (et);
- if (fd17 < 0)
- G_fatal_error ("unable to create raster map %s", et);
- }
- if (tc != NULL)
- {
- cell13 = G_allocate_f_raster_buf();
- /* if (ts == 1) {
- sprintf(buf,"%s.%.*d",tc,ndigit,tt);
- tc0 = G_store(buf);
- fd13 = G_open_fp_cell_new (tc0);
- }
- else*/
- fd13 = G_open_fp_cell_new (tc);
+ if (et != NULL) {
+ cell17 = G_allocate_f_raster_buf();
+ /* if (ts == 1) {
+ sprintf(buf,"%s.%.*d",et,ndigit,tt);
+ et0 = G_store(buf);
+ fd17 = G_open_fp_cell_new (et0);
+ }
+ else */
+ fd17 = G_open_fp_cell_new(et);
- if (fd13 < 0)
- G_fatal_error ("unable to create raster map %s", tc);
- }
+ if (fd17 < 0)
+ G_fatal_error("unable to create raster map %s", et);
+ }
+ if (tc != NULL) {
+ cell13 = G_allocate_f_raster_buf();
+ /* if (ts == 1) {
+ sprintf(buf,"%s.%.*d",tc,ndigit,tt);
+ tc0 = G_store(buf);
+ fd13 = G_open_fp_cell_new (tc0);
+ }
+ else */
+ fd13 = G_open_fp_cell_new(tc);
- if(G_set_window (&cellhd) < 0)
- G_fatal_error ("G_set_window");
+ if (fd13 < 0)
+ G_fatal_error("unable to create raster map %s", tc);
+ }
- if (my != G_window_rows())
- G_fatal_error ("OOPS: rows changed from %d to %d\n", mx, G_window_rows());
- if (mx != G_window_cols())
- G_fatal_error ("OOPS: cols changed from %d to %d\n", my, G_window_cols());
- for(iarc=0;iarc<my;iarc++)
- {
- i=my-iarc-1;
- if(et!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || er[i][j] == UNDEF)
- G_set_f_null_value (cell17+j,1);
- else {
- cell17[j]=(FCELL) er[i][j];/* add conv? */
- etmax = amax1(etmax, er[i][j]);
- etmin = amin1(etmin, er[i][j]);
- }
- }
- G_put_f_raster_row (fd17, cell17);
- }
+ if (G_set_window(&cellhd) < 0)
+ G_fatal_error("G_set_window");
- if(tc!=NULL)
- {
- for(j=0;j<mx;j++)
- {
- if(zz[i][j] == UNDEF || sigma[i][j] == UNDEF || si[i][j] == UNDEF)
- G_set_f_null_value (cell13+j,1);
+ if (my != G_window_rows())
+ G_fatal_error("OOPS: rows changed from %d to %d\n", mx,
+ G_window_rows());
+ if (mx != G_window_cols())
+ G_fatal_error("OOPS: cols changed from %d to %d\n", my,
+ G_window_cols());
+
+ for (iarc = 0; iarc < my; iarc++) {
+ i = my - iarc - 1;
+ if (et != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || er[i][j] == UNDEF)
+ G_set_f_null_value(cell17 + j, 1);
else {
- if (sigma[i][j] == 0.) trc = 0.;
- else
+ cell17[j] = (FCELL) er[i][j]; /* add conv? */
+ etmax = amax1(etmax, er[i][j]);
+ etmin = amin1(etmin, er[i][j]);
+ }
+ }
+ G_put_f_raster_row(fd17, cell17);
+ }
+
+ if (tc != NULL) {
+ for (j = 0; j < mx; j++) {
+ if (zz[i][j] == UNDEF || sigma[i][j] == UNDEF ||
+ si[i][j] == UNDEF)
+ G_set_f_null_value(cell13 + j, 1);
+ else {
+ if (sigma[i][j] == 0.)
+ trc = 0.;
+ else
trc = si[i][j] / sigma[i][j];
- cell13[j]=(FCELL) trc;
- /* gsmax = amax1(gsmax, trc); */
+ cell13[j] = (FCELL) trc;
+ /* gsmax = amax1(gsmax, trc); */
}
- }
- G_put_f_raster_row (fd13, cell13);
- }
+ }
+ G_put_f_raster_row(fd13, cell13);
}
+ }
- if(tc!=NULL)
- G_close_cell (fd13);
+ if (tc != NULL)
+ G_close_cell(fd13);
- if(et!=NULL)
- G_close_cell (fd17);
+ if (et != NULL)
+ G_close_cell(fd17);
- if (et!=NULL) {
+ if (et != NULL) {
- G_init_colors(&colors);
+ G_init_colors(&colors);
- dat1 = (FCELL) etmax;
- dat2 = (FCELL) 0.1;
- G_add_f_raster_color_rule(&dat1,0,0,0,&dat2,0,0,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.01;
- G_add_f_raster_color_rule(&dat1,0,0,255,&dat2,0,191,191,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.0001;
- G_add_f_raster_color_rule(&dat1,0,191,191,&dat2,170,255,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) 0.;
- G_add_f_raster_color_rule(&dat1,170,255,255,&dat2,255,255,255,&colors);
- dat1 = dat2;
- dat2 = (FCELL) -0.0001;
- G_add_f_raster_color_rule(&dat1,255,255,255,&dat2,255,255,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) -0.01;
- G_add_f_raster_color_rule(&dat1,255,255,0,&dat2,255,127,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) -0.1;
- G_add_f_raster_color_rule(&dat1,255,127,0,&dat2,255,0,0,&colors);
- dat1 = dat2;
- dat2 = (FCELL) etmin;
- G_add_f_raster_color_rule(&dat1,255,0,0,&dat2,255,0,255,&colors);
+ dat1 = (FCELL) etmax;
+ dat2 = (FCELL) 0.1;
+ G_add_f_raster_color_rule(&dat1, 0, 0, 0, &dat2, 0, 0, 255, &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.01;
+ G_add_f_raster_color_rule(&dat1, 0, 0, 255, &dat2, 0, 191, 191,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.0001;
+ G_add_f_raster_color_rule(&dat1, 0, 191, 191, &dat2, 170, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) 0.;
+ G_add_f_raster_color_rule(&dat1, 170, 255, 255, &dat2, 255, 255, 255,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.0001;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.01;
+ G_add_f_raster_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) - 0.1;
+ G_add_f_raster_color_rule(&dat1, 255, 127, 0, &dat2, 255, 0, 0,
+ &colors);
+ dat1 = dat2;
+ dat2 = (FCELL) etmin;
+ G_add_f_raster_color_rule(&dat1, 255, 0, 0, &dat2, 255, 0, 255,
+ &colors);
-/* if (ts == 1) {
- if ((mapst = G_find_file("cell", et0, "")) == NULL)
- G_fatal_error ("cannot find file %s", et0);
- G_write_colors(et0, mapst, &colors);
- G_quantize_fp_map_range(et0,mapst,(FCELL)etmin,(FCELL)etmax,(CELL)etmin,(CELL)etmax);
- G_free_colors(&colors);
+ /* if (ts == 1) {
+ if ((mapst = G_find_file("cell", et0, "")) == NULL)
+ G_fatal_error ("cannot find file %s", et0);
+ G_write_colors(et0, mapst, &colors);
+ G_quantize_fp_map_range(et0,mapst,(FCELL)etmin,(FCELL)etmax,(CELL)etmin,(CELL)etmax);
+ G_free_colors(&colors);
+ }
+ else { */
+
+ if ((mapst = G_find_file("cell", et, "")) == NULL)
+ G_fatal_error("cannot find file %s", et);
+ G_write_colors(et, mapst, &colors);
+ G_quantize_fp_map_range(et, mapst, (FCELL) etmin, (FCELL) etmax,
+ (CELL) etmin, (CELL) etmax);
+ G_free_colors(&colors);
+ /* } */
}
- else { */
- if ((mapst = G_find_file("cell", et, "")) == NULL)
- G_fatal_error ("cannot find file %s", et);
- G_write_colors(et, mapst, &colors);
- G_quantize_fp_map_range(et, mapst,(FCELL)etmin,(FCELL)etmax,(CELL)etmin,(CELL)etmax);
- G_free_colors(&colors);
- /* } */
- }
-
- return 1;
+ return 1;
}
-
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/random.c
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/random.c 2008-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/random.c 2008-09-13 08:32:00 UTC (rev 33426)
@@ -15,7 +15,7 @@
seed.is1 = irand1;
seed.is2 = irand2;
return 0;
-}
+}
int seedg(long int irand1, long int irand2)
@@ -23,7 +23,7 @@
irand1 = seed.is1;
irand2 = seed.is2;
return 0;
-}
+}
double ulec(void)
@@ -35,15 +35,15 @@
long int k, iz;
-/* uniform random number generator (combined type) */
-/* P. L'Ecuyer, Commun. ACM, 31(1988)742 */
-/* portable (32 bits arithmetics) */
+ /* uniform random number generator (combined type) */
+ /* P. L'Ecuyer, Commun. ACM, 31(1988)742 */
+ /* portable (32 bits arithmetics) */
k = seed.is1 / 53668;
seed.is1 -= k * 53668;
seed.is1 = seed.is1 * 40014 - k * 12211;
-/* is1=40014*(is1-k*53668)-k*12211 */
+ /* is1=40014*(is1-k*53668)-k*12211 */
if (seed.is1 < 0) {
seed.is1 += 2147483563;
}
@@ -51,7 +51,7 @@
k = seed.is2 / 52774;
seed.is2 -= k * 52774;
seed.is2 = seed.is2 * 40692 - k * 3791;
-/* is2=40692*(is2-k*52774)-k*3791 */
+ /* is2=40692*(is2-k*52774)-k*3791 */
if (seed.is2 < 0) {
seed.is2 += 2147483399;
}
@@ -60,10 +60,10 @@
if (iz < 0) {
iz += 2147483562;
}
- ret_val = (double) iz * 4.656613e-10;
+ ret_val = (double)iz *4.656613e-10;
return ret_val;
-} /* ulec */
+} /* ulec */
double gasdev(void)
@@ -77,22 +77,22 @@
double ret_val;
/* Local variables */
- double r=0., vv1, vv2, fac;
+ double r = 0., vv1, vv2, fac;
if (iset == 0) {
while (r >= 1. || r == 0.) {
- vv1 = ulec() * 2. - 1.;
- vv2 = ulec() * 2. - 1.;
- r = vv1 * vv1 + vv2 * vv2;
+ vv1 = ulec() * 2. - 1.;
+ vv2 = ulec() * 2. - 1.;
+ r = vv1 * vv1 + vv2 * vv2;
}
fac = sqrt(log(r) * -2. / r);
gset = vv1 * fac;
ret_val = vv2 * fac;
iset = 1;
- } else {
+ }
+ else {
ret_val = gset;
iset = 0;
}
return ret_val;
-} /* gasdev */
-
+} /* gasdev */
Modified: grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h
===================================================================
--- grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h 2008-09-13 08:28:10 UTC (rev 33425)
+++ grass/branches/releasebranch_6_3/raster/simwe/simlib/waterglobs.h 2008-09-13 08:32:00 UTC (rev 33426)
@@ -11,9 +11,11 @@
#define GLOBAL extern
#endif
-GLOBAL FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps, *fdmanin, *fddepth, *fddisch, *fderr, *fdoutwalk,*fdwalkers;
-GLOBAL FILE *fdwdepth,*fddetin,*fdtranin,*fdtauin, *fdtc, *fdet, *fdconc, *fdflux,*fderdep;
-GLOBAL FILE *fdsfile,*fw;
+GLOBAL FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps,
+ *fdmanin, *fddepth, *fddisch, *fderr, *fdoutwalk, *fdwalkers;
+GLOBAL FILE *fdwdepth, *fddetin, *fdtranin, *fdtauin, *fdtc, *fdet, *fdconc,
+ *fdflux, *fderdep;
+GLOBAL FILE *fdsfile, *fw;
GLOBAL char *elevin;
GLOBAL char *dxin;
@@ -47,9 +49,11 @@
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,*rainval,*maninval,*infilval;
+ 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, *rainval, *maninval,
+ *infilval;
} parm;
@@ -78,53 +82,53 @@
GLOBAL int npoints_alloc;
GLOBAL int input_data(void);
-GLOBAL int seeds(long int,long int);
-GLOBAL int seedg(long int,long int);
+GLOBAL int seeds(long int, long int);
+GLOBAL int seedg(long int, long int);
GLOBAL int grad_check(void);
GLOBAL void erod(double **);
GLOBAL void main_loop(void);
-GLOBAL int output_data(int,double);
+GLOBAL int output_data(int, double);
GLOBAL int output_et(void);
GLOBAL double ulec(void);
GLOBAL double gasdev(void);
-GLOBAL double amax1(double,double);
-GLOBAL double amin1(double,double);
-GLOBAL int min(int,int);
-GLOBAL int max(int,int);
+GLOBAL double amax1(double, double);
+GLOBAL double amin1(double, double);
+GLOBAL int min(int, int);
+GLOBAL int max(int, int);
GLOBAL double xmin, ymin, xmax, ymax;
GLOBAL double mayy, miyy, maxx, mixx;
GLOBAL int mx, my;
GLOBAL int mx2, my2;
-GLOBAL double bxmi,bymi,bxma,byma,bresx,bresy;
+GLOBAL double bxmi, bymi, bxma, byma, bresx, bresy;
GLOBAL int maxwab;
-GLOBAL double step,conv;
+GLOBAL double step, conv;
GLOBAL double frac;
GLOBAL double bxmi, bymi;
GLOBAL float **zz, **cchez;
GLOBAL double **v1, **v2, **slope;
-GLOBAL double **gama, **gammas,**si,**inf,**sigma;
-GLOBAL float **dc,**tau,**er, **ct, **trap;
-GLOBAL float **dif;
+GLOBAL double **gama, **gammas, **si, **inf, **sigma;
+GLOBAL float **dc, **tau, **er, **ct, **trap;
+GLOBAL float **dif;
-GLOBAL double vavg[MAXW][2], stack[MAXW][3],w[MAXW][3];
+GLOBAL double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
GLOBAL int iflag[MAXW];
GLOBAL double hbeta;
GLOBAL int ldemo;
-GLOBAL double hhmax, sisum,vmean;
-GLOBAL double infsum,infmean;
+GLOBAL double hhmax, sisum, vmean;
+GLOBAL double infsum, infmean;
GLOBAL int maxw, maxwa, nwalk;
GLOBAL double rwalk, bresx, bresy, xrand, yrand;
GLOBAL double stepx, stepy, xp0, yp0;
-GLOBAL double chmean, si0, deltap, deldif, cch, hhc,halpha;
+GLOBAL double chmean, si0, deltap, deldif, cch, hhc, halpha;
GLOBAL double eps;
GLOBAL int maxwab, nstack;
GLOBAL int iterout, mx2o, my2o;
-GLOBAL int miter,nwalka,lwwfin;
+GLOBAL int miter, nwalka, lwwfin;
GLOBAL double timec;
GLOBAL int ts, timesec;
More information about the grass-commit
mailing list