MFD stable, several enhancements

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@35412 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.watershed/front/main.c

@@ -20,7 +20,7 @@
 #include <grass/gis.h>
 #include <grass/glocale.h>
 
-int write_hist(char *, char *, char *, int);
+int write_hist(char *, char *, char *, int, int);
 
 int main(int argc, char *argv[])
 {
@@ -62,7 +62,8 @@ int main(int argc, char *argv[])
 
     opt2 = G_define_option();
     opt2->key = "depression";
-    opt2->description = _("Input map: locations of real depressions");
+    opt2->label = _("Input map: locations of real depressions");
+    opt2->description = _("All non-NULL cells are considered as real depressions");
     opt2->required = NO;
     opt2->type = TYPE_STRING;
     opt2->gisprompt = "old,cell,raster";
@@ -87,8 +88,10 @@ int main(int argc, char *argv[])
 
     opt5 = G_define_option();
     opt5->key = "blocking";
-    opt5->description =
+    opt5->label =
 	_("Input map: terrain blocking overland surface flow, for USLE");
+    opt5->description =
+	_("All non-NULL cells are considered as blocking terrain");
     opt5->required = NO;
     opt5->type = TYPE_STRING;
     opt5->gisprompt = "old,cell,raster";
@@ -232,18 +235,20 @@ int main(int argc, char *argv[])
     }
 
     err = 0;
-    /* basin and basin.thresh */
+    /* basin and basin threshold */
     err += (opt10->answer != NULL && opt6->answer == NULL);
-    /* stream and basin.thresh */
+    /* stream and basin threshold */
     err += (opt11->answer != NULL && opt6->answer == NULL);
-    /* half.basin and basin.thresh */
+    /* half_basin and basin threshold */
     err += (opt12->answer != NULL && opt6->answer == NULL);
-    /* slope and basin.thresh */
+    /* LS factor and basin threshold */
+    err += (opt14->answer != NULL && opt6->answer == NULL);
+    /* S factor and basin threshold */
     err += (opt15->answer != NULL && opt6->answer == NULL);
 
     if (err) {
 	G_message(_("Sorry, if any of the following options are set:\n"
-		    "    basin, stream, half.basin, slope, or lS\n"
+		    "    basin, stream, half_basin, length_slope, or slope_steepness\n"
 		    "    you MUST provide a value for the basin "
 		    "threshold parameter."));
 	G_usage();
@@ -388,39 +393,41 @@ int main(int argc, char *argv[])
     if (opt8->answer)
 	write_hist(opt8->answer,
 		   "Watershed accumulation: overland flow that traverses each cell",
-		   opt1->answer, flag_seg->answer);
+		   opt1->answer, flag_seg->answer, flag_sfd->answer);
     if (opt9->answer)
 	write_hist(opt9->answer,
 		   "Watershed drainage direction (divided by 45deg)",
-		   opt1->answer, flag_seg->answer);
+		   opt1->answer, flag_seg->answer, flag_sfd->answer);
     if (opt10->answer)
 	write_hist(opt10->answer,
-		   "Watershed basins", opt1->answer, flag_seg->answer);
+		   "Watershed basins", opt1->answer, flag_seg->answer, 
+		   flag_sfd->answer);
     if (opt11->answer)
 	write_hist(opt11->answer,
 		   "Watershed stream segments", opt1->answer,
-		   flag_seg->answer);
+		   flag_seg->answer, flag_sfd->answer);
     if (opt12->answer)
 	write_hist(opt12->answer,
-		   "Watershed half-basins", opt1->answer, flag_seg->answer);
+		   "Watershed half-basins", opt1->answer, flag_seg->answer, 
+		   flag_sfd->answer);
     if (opt13->answer)
 	write_hist(opt13->answer,
 		   "Watershed visualization map (filtered accumulation map)",
-		   opt1->answer, flag_seg->answer);
+		   opt1->answer, flag_seg->answer, flag_sfd->answer);
     if (opt14->answer)
 	write_hist(opt14->answer,
 		   "Watershed slope length and steepness (LS) factor",
-		   opt1->answer, flag_seg->answer);
+		   opt1->answer, flag_seg->answer, flag_sfd->answer);
     if (opt15->answer)
 	write_hist(opt15->answer,
 		   "Watershed slope steepness (S) factor",
-		   opt1->answer, flag_seg->answer);
+		   opt1->answer, flag_seg->answer, flag_sfd->answer);
 
     exit(ret);
 }
 
 /* record map history info */
-int write_hist(char *map_name, char *title, char *source_name, int mode)
+int write_hist(char *map_name, char *title, char *source_name, int mode, int sfd)
 {
     struct History history;
 
@@ -430,8 +437,10 @@ int write_hist(char *map_name, char *title, char *source_name, int mode)
     strncpy(history.datsrc_1, source_name, RECORD_LEN);
     history.datsrc_1[RECORD_LEN - 1] = '\0';	/* strncpy() doesn't null terminate if maxfill */
     sprintf(history.edhist[0],
-	    "Processing mode: %s", mode ? "Segmented" : "All in RAM");
-    history.edlinecnt = 1;
+	    "Processing mode: %s", sfd ? "SFD (D8)" : "MFD");
+    sprintf(history.edhist[1],
+	    "Memory mode: %s", mode ? "Segmented" : "All in RAM");
+    history.edlinecnt = 2;
     G_command_history(&history);
 
     return G_write_history(map_name, &history);

raster/r.watershed/front/r.watershed.html

@@ -1,7 +1,7 @@
 <h2>DESCRIPTION</h2>
 
 <em>r.watershed</em> generates a set of maps indicating:
-1) the location of watershed basins, and
+1) flow accumulation, drainage direction, the location of streams and watershed basins, and
 2) the LS and S factors of the Revised Universal Soil Loss Equation (RUSLE).
 
 <p>
@@ -24,13 +24,13 @@ lumped-parameter hydrologic model program.  The non-interactive version of
 
 <dd>Without this flag set, the entire analysis is run in memory
 maintained by the operating system.  This can be limiting, but is
-relatively fast.  Setting the flag causes the program to manage memory
+very fast.  Setting the flag causes the program to manage memory
 on disk which allows larger maps to be processed but is considerably
 slower.
 
 <dt><em>-s</em> 
 
-<dd>Use single flow direction (SFD) instead of multiple flow direction (MFD).
+<dd>Use single flow direction (SFD, D8) instead of multiple flow direction (MFD).
 MFD is enabled by default.
 
 <dt><em>-4</em> 
@@ -59,7 +59,8 @@ flow is less widely distributed, becoming more similar to SFD.
 
 <dd>Input map:  Map layer of actual depressions or sinkholes in the
 landscape that are large enough to slow and store surface runoff from 
-a storm event.  Any non-zero values indicate depressions.
+a storm event.  Any non-NULL (not nodata) values indicate depressions. 
+Water will flow into depressions, but not out of depressions.
 
 <dt><em>flow</em> 
 
@@ -67,7 +68,7 @@ a storm event.  Any non-zero values indicate depressions.
 amount of overland flow units that each cell will contribute to the
 watershed basin model.  Overland flow units represent the amount of
 overland flow each cell contributes to surface flow.  If omitted, a
-value of one (1) is assumed. The algorithm is D8 flow accumulation.
+value of one (1) is assumed.
 
 <dt><em>disturbed_land</em> 
 
@@ -80,7 +81,7 @@ assumes no disturbed land.  This input is used for the RUSLE calculations.
 
 <dd>Input map: terrain that will block overland surface flow.  Terrain
 that will block overland surface flow and restart the slope length
-for the RUSLE.  Any non-zero values indicate blocking terrain.
+for the RUSLE.  Any non-NULL (not nodata) values indicate blocking terrain.
 
 <dt><em>threshold</em> 
 
@@ -116,8 +117,8 @@ calculated accurately.
 <dd>Output map: drainage direction.  Provides the "aspect" for each
 cell.  Multiplying positive values by 45 will give the direction in
 degrees that the surface runoff will travel from that cell.  The
-value -1 indicates that the cell is a depression area (defined by
-the depression input map).  Other negative values indicate that
+value 0 (zero) indicates that the cell is a depression area (defined by
+the depression input map).  Negative values indicate that
 surface runoff is leaving the boundaries of the current geographic
 region.  The absolute value of these negative cells indicates the
 direction of flow.
@@ -154,22 +155,19 @@ are given the value of the <em>threshold</em> parameter.
 
 <dt><em>length_slope</em> 
 
-<dd>Output map: slope length and steepness (LS) factor.  Contains the LS
-factor for the Revised Universal Soil Loss Equation.  Equations taken
-from <em>Revised Universal Soil Loss Equation for Western Rangelands</em>
-(Weltz et al. 1987).
-Since the LS factor is a small number, it is multiplied by 100 for the
-GRASS output map.
+<dd>Output map: slope length and steepness (LS) factor for the Revised 
+Universal Soil Loss Equation (RUSLE).  Equations taken from <em>Revised 
+Universal Soil Loss Equation for Western Rangelands</em>
+(Weltz et al. 1987). Since the LS factor is a small number (usually less 
+than one), the GRASS output map is of type DCELL.
 
 <dt><em>slope_steepness</em> 
 
-<dd>Output map: slope steepness (S) factor for RUSLE.
-Contains the revised S factor for the Universal Soil
-Loss Equation.  Equations taken from article entitled
+<dd>Output map: slope steepness (S) factor for the Universal Soil
+Loss Equation (RUSLE).  Equations taken from article entitled
 <em>Revised Slope Steepness Factor for the Universal Soil
-Loss Equation</em> (McCool et al. 1987).  Since the S factor
-is a small number (usually less than one), it is multiplied
-by 100 for the GRASS output map layer.
+Loss Equation</em> (McCool et al. 1987).  Since the LS factor is a small 
+number (usually less than one), the GRASS output map is of type DCELL.
 </dd>
 </dl>
 
@@ -208,13 +206,13 @@ single flow direction (SFD, D8) and multiple flow direction (MFD). With
 MFD, water flow is distributed to all neighbouring cells with lower 
 elevation, using slope towards neighbouring cells as a weighing factor 
 for proportional distribution. The A<sup>T</sup> least-cost path is 
-always included and assigned the maxmimum observed weighing factor. 
-As a result, depressions and obstacles are overflown with a gracefull 
-flow convergence before the overflow. The convergence factor causes flow 
-accumulation to converge more strongly with higher values. The supported 
-range is 1 to 10, recommended is a convergence factor of 5. If many 
-small sliver basins are created with MFD, setting the convergence factor 
-to a higher value can reduce the amount of small sliver basins.
+always included. As a result, depressions and obstacles are overflown 
+with a gracefull flow convergence before the overflow. The convergence 
+factor causes flow accumulation to converge more strongly with higher 
+values. The supported range is 1 to 10, recommended is a convergence 
+factor of 5 (Holmgren, 1994). If many small sliver basins are created 
+with MFD, setting the convergence factor to a higher value can reduce 
+the amount of small sliver basins.
 <br>See example below with the North Carolina dataset for using MFD mode.
 
 <h4>In-memory mode and disk swap mode</h4>
@@ -222,12 +220,17 @@ There are two versions of this program: <em>ram</em> and <em>seg</em>.
 <em>ram</em> is used by default, <em>seg</em> can be used by setting 
 the <em>-m</em> flag.
 <br>
+The <em>ram</em> version requires a maximum of 31 MB of RAM for 1 million 
+cells. Together with the amount of system memory (RAM) available, this 
+value can be used to estimate whether the current region can be 
+processed with the <em>ram</em> version.
+<br>
 The <em>ram</em> version uses virtual memory managed by the operating
 system to store all the data structures and is faster than the <em>seg</em>
 version; <em>seg</em> uses the GRASS segmentation library which manages 
-data in disk files. Thus <em>seg</em> uses much less system memory (RAM) 
-allowing other processes to operate on the same CPU, even when the 
-current geographic region is huge.
+data in disk files. <em>seg</em> uses only as much system memory (RAM) as 
+specified with the <em>memory</em> option, allowing other processes to 
+operate on the same system, even when the current geographic region is huge.
 <br>
 Due to memory requirements of both programs, it is quite easy to run out of
 memory when working with huge map regions. If the <em>ram</em> version runs
@@ -280,15 +283,16 @@ tributaries flowing into it.
 <p>
 The <em>r.watershed</em> program does not require the user to have the
 current geographic region filled with elevation values.  Areas without
-elevation data MUST be masked out, by creating a raster map (or raster
-reclassification) named <tt>MASK</tt>.  Areas
-masked out will be treated as if they are off the edge of the region.
-MASKs will reduce the memory necessary to run the program.  Masking out 
-unimportant areas can significantly reduce processing time if the watersheds 
-of interest occupy a small percentage of the overall area.
+elevation data (masked or NULL cells) are ignored. It is NOT necessary to
+create a raster map (or raster reclassification) named <tt>MASK</tt> for 
+NULL cells.  Areas without elevation data will be treated as if they are 
+off the edge of the region. Such areas will reduce the memory necessary 
+to run the program.  Masking out unimportant areas can significantly 
+reduce processing time if the watersheds of interest occupy a small 
+percentage of the overall area.
 
 <p>
-Zero data values will be treated as elevation data (not no_data).
+Zero (0) and negative data values will be treated as elevation data (not no_data).
 
 <h4>Further processing of output layers</h4>
 <p>

raster/r.watershed/ram/Gwater.h

@@ -34,7 +34,7 @@
 
 #define POINT       struct points
 POINT {
-    SHORT r, c, downr, downc;
+    SHORT r, c; /* , downr, downc */
     int nxt;
 };
 
@@ -87,7 +87,7 @@ CELL def_basin(int, int, CELL, double, CELL);
 
 /* do_astar.c */
 int do_astar(void);
-int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int add_pt(SHORT, SHORT, CELL, CELL);
 int drop_pt(void);
 int sift_up(int, CELL);
 double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);

raster/r.watershed/ram/close_maps.c

@@ -6,21 +6,24 @@
 
 int close_maps(void)
 {
-    struct Colors colors, logcolors;
+    struct Colors colors;
     int value, r, c, fd;
     CELL *buf = NULL;
     DCELL *dbuf = NULL;
     struct FPRange accRange;
     DCELL min, max;
     DCELL clr_min, clr_max;
+    DCELL sum, sum_sqr, stddev, lstddev, dvalue;
 
-    if (wat_flag || asp_flag || dis_flag || ls_flag || sl_flag || sg_flag)
+    if (asp_flag || dis_flag)
 	buf = G_allocate_cell_buf();
-    dbuf = G_allocate_d_raster_buf();
+    if (wat_flag || ls_flag || sl_flag || sg_flag)
+	dbuf = G_allocate_d_raster_buf();
     G_free(alt);
     if (ls_flag || sg_flag)
 	G_free(r_h);
 
+    sum = sum_sqr = stddev = 0.0;
     if (wat_flag) {
 	fd = G_open_raster_new(wat_name, DCELL_TYPE);
 	if (fd < 0) {
@@ -28,65 +31,92 @@ int close_maps(void)
 	}
 	else {
 	    for (r = 0; r < nrows; r++) {
+		G_set_d_null_value(dbuf, ncols);	/* reset row to all NULL */
 		for (c = 0; c < ncols; c++) {
-		    dbuf[c] = wat[SEG_INDEX(wat_seg, r, c)];
+		    dvalue = wat[SEG_INDEX(wat_seg, r, c)];
+		    if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+			dbuf[c] = dvalue;
+			dvalue = ABS(dvalue);
+			sum += dvalue;
+			sum_sqr += dvalue * dvalue;
+		    }
 		}
 		G_put_raster_row(fd, dbuf, DCELL_TYPE);
 	    }
 	    if (G_close_cell(fd) < 0)
 		G_warning(_("Close failed."));
 
+	    stddev =
+		sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
+	    G_debug(1, "stddev: %f", stddev);
+
 	    /* set nice color rules: yellow, green, cyan, blue, black */
+
+	    lstddev = log(stddev);
+
 	    G_read_fp_range(wat_name, this_mapset, &accRange);
 	    min = max = 0;
 	    G_get_fp_range_min_max(&accRange, &min, &max);
 
 	    G_init_colors(&colors);
-	    if (ABS(min) > max) {
-		clr_min = 1;
-		clr_max = 0.25 * max;
-		G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
-					  255, 0, &colors);
-		clr_min = 0.25 * max;
-		clr_max = 0.5 * max;
-		G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
-					  255, 255, &colors);
-		clr_min = 0.5 * max;
-		clr_max = 0.75 * max;
-		G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+
+	    if (min < 0) {
+		if (min < (-stddev - 1)) {
+		    clr_min = min;
+		    clr_max = -stddev - 1;
+		    G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+					      0, 0, &colors);
+		}
+		clr_min = -stddev - 1.;
+		clr_max = -1. * exp(lstddev * 0.75);
+		G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
 					  0, 255, &colors);
-		clr_min = 0.75 * max;
-		clr_max = max;
-		G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
-					  0, &colors);
-		max = -max;
-	    }
-	    else {
-		min = ABS(min);
-		clr_min = 1;
-		clr_max = 0.25 * min;
-		G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
-					  255, 0, &colors);
-		clr_min = 0.25 * min;
-		clr_max = 0.5 * min;
-		G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
+		clr_min = clr_max;
+		clr_max = -1. * exp(lstddev * 0.5);
+		G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0,
 					  255, 255, &colors);
-		clr_min = 0.5 * min;
-		clr_max = 0.75 * min;
+		clr_min = clr_max;
+		clr_max = -1. * exp(lstddev * 0.35);
 		G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
-					  0, 255, &colors);
-		clr_min = 0.75 * min;
-		clr_max = min;
-		G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
+					  255, 0, &colors);
+		clr_min = clr_max;
+		clr_max = -1.;
+		G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 255,
+					  255, 0, &colors);
+	    }
+	    clr_min = -1.;
+	    clr_max = 1.;
+	    G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 255,
+				      255, 0, &colors);
+	    clr_min = 1;
+	    clr_max = exp(lstddev * 0.35);
+	    G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
+				      255, 0, &colors);
+	    clr_min = clr_max;
+	    clr_max = exp(lstddev * 0.5);
+	    G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
+				      255, 255, &colors);
+	    clr_min = clr_max;
+	    clr_max = exp(lstddev * 0.75);
+	    G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+				      0, 255, &colors);
+	    clr_min = clr_max;
+	    clr_max = stddev + 1.;
+	    G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
+				      0, &colors);
+
+	    if (max > 0 && max > stddev + 1) {
+		clr_min = stddev + 1;
+		clr_max = max;
+		G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0, 0,
 					  0, &colors);
 	    }
-	    G_abs_log_colors(&logcolors, &colors, 5);
-	    G_add_d_raster_color_rule(&min, 0, 0, 0, &max, 0, 0, 0,
-				      &logcolors);
-	    G_write_colors(wat_name, this_mapset, &logcolors);
+	    G_write_colors(wat_name, this_mapset, &colors);
 	}
     }
 
+    /* TODO: elevation == NULL -> drainage direction == NULL (wat == 0 where ele == NULL) */
+    /* keep drainage direction == 0 for real depressions */
     if (asp_flag) {
 	fd = G_open_cell_new(asp_name);
 	if (fd < 0) {
@@ -94,6 +124,7 @@ int close_maps(void)
 	}
 	else {
 	    for (r = 0; r < nrows; r++) {
+		G_set_c_null_value(buf, ncols);	/* reset row to all NULL */
 		for (c = 0; c < ncols; c++) {
 		    buf[c] = asp[SEG_INDEX(asp_seg, r, c)];
 		}
@@ -146,16 +177,16 @@ int close_maps(void)
     G_free(wat);
 
     if (ls_flag) {
-	fd = G_open_cell_new(ls_name);
+	fd = G_open_raster_new(ls_name, DCELL_TYPE);
 	if (fd < 0) {
-	    G_warning(_("unable to open new L map layer."));
+	    G_warning(_("unable to open new LS factor map layer."));
 	}
 	else {
 	    for (r = 0; r < nrows; r++) {
 		for (c = 0; c < ncols; c++) {
-		    buf[c] = l_s[SEG_INDEX(l_s_seg, r, c)] + .5;
+		    dbuf[c] = l_s[SEG_INDEX(l_s_seg, r, c)];
 		}
-		G_put_raster_row(fd, buf, CELL_TYPE);
+		G_put_raster_row(fd, dbuf, DCELL_TYPE);
 	    }
 	    if (G_close_cell(fd) < 0)
 		G_warning(_("Close failed."));
@@ -164,18 +195,18 @@ int close_maps(void)
     }
 
     if (sl_flag) {
-	fd = G_open_cell_new(sl_name);
+	fd = G_open_raster_new(sl_name, DCELL_TYPE);
 	if (fd < 0) {
 	    G_warning(_("unable to open new slope length map layer."));
 	}
 	else {
 	    for (r = 0; r < nrows; r++) {
 		for (c = 0; c < ncols; c++) {
-		    buf[c] = s_l[SEG_INDEX(s_l_seg, r, c)] + .5;
-		    if (buf[c] > max_length)
-			buf[c] = max_length + .5;
+		    dbuf[c] = s_l[SEG_INDEX(s_l_seg, r, c)];
+		    if (dbuf[c] > max_length)
+			dbuf[c] = max_length;
 		}
-		G_put_raster_row(fd, buf, CELL_TYPE);
+		G_put_raster_row(fd, dbuf, DCELL_TYPE);
 	    }
 	    if (G_close_cell(fd) < 0)
 		G_warning(_("Close failed."));
@@ -186,16 +217,16 @@ int close_maps(void)
 	G_free(s_l);
 
     if (sg_flag) {
-	fd = G_open_cell_new(sg_name);
+	fd = G_open_raster_new(sg_name, DCELL_TYPE);
 	if (fd < 0) {
-	    G_warning(_("unable to open new S map layer."));
+	    G_warning(_("unable to open new S factor map layer."));
 	}
 	else {
 	    for (r = 0; r < nrows; r++) {
 		for (c = 0; c < ncols; c++) {
-		    buf[c] = s_g[SEG_INDEX(s_g_seg, r, c)] * 100 + .5;
+		    dbuf[c] = s_g[SEG_INDEX(s_g_seg, r, c)];
 		}
-		G_put_raster_row(fd, buf, CELL_TYPE);
+		G_put_raster_row(fd, dbuf, DCELL_TYPE);
 	    }
 	    if (G_close_cell(fd) < 0)
 		G_warning(_("Close failed."));

raster/r.watershed/ram/close_maps2.c

@@ -112,7 +112,7 @@ int close_array_seg(void)
 	G_write_colors(bas_name, this_mapset, &colors);
     }
 
-    /* half.basins map */
+    /* half_basins map */
     if (haf_flag) {
 	map_fd = G_open_cell_new(haf_name);
 	for (r = 0; r < nrows; r++) {

raster/r.watershed/ram/do_astar.c

@@ -36,7 +36,7 @@ int do_astar(void)
 	/* equivalent to first_astar = point->next in old code */
 	drop_pt();
 
-	/* can go, dragged on from old code */
+	/* can go, dragged on from old code: replace first_astar with heap_index[1] in line 22 */
 	first_astar = heap_index[1];
 
 	/* downhill path for flow accumulation is set here */
@@ -51,7 +51,6 @@ int do_astar(void)
 
 	index_doer = SEG_INDEX(alt_seg, r, c);
 	alt_val = alt[index_doer];
-	wat_val = wat[index_doer];
 
 	FLAG_SET(worked, r, c);
 
@@ -69,24 +68,24 @@ int do_astar(void)
 		if (in_val == 0) {
 		    alt_up = alt[index_up];
 		    /* flow direction is set here */
-		    add_pt(upr, upc, r, c, alt_up, alt_val);
+		    add_pt(upr, upc, alt_up, alt_val);
 		    drain_val = drain[upr - r + 1][upc - c + 1];
 		    asp[index_up] = drain_val;
-
 		}
 		else {
-		    /* check if neighbour has not been worked on,
-		     * update values for asp and wat */
+		    /* check if neighbour has been worked on,
+		     * if not, update values for asp and wat */
 		    in_val = FLAG_GET(worked, upr, upc);
 		    if (in_val == 0) {
 			asp_up = asp[index_up];
-			if (asp_up < -1) {
+			if (asp_up < 0) {
 			    asp[index_up] = drain[upr - r + 1][upc - c + 1];
 
+			    wat_val = wat[index_doer];
 			    if (wat_val > 0)
 				wat[index_doer] = -wat_val;
 
-			    replace(upr, upc, r, c);	/* alt_up used to be */
+			    /* replace(upr, upc, r, c); */	/* alt_up used to be */
 			}
 		    }
 		}
@@ -104,7 +103,7 @@ int do_astar(void)
 }
 
 /* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
 {
 
     FLAG_SET(in_list, r, c);
@@ -120,8 +119,8 @@ int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 
     astar_pts[nxt_avail_pt].r = r;
     astar_pts[nxt_avail_pt].c = c;
-    astar_pts[nxt_avail_pt].downr = downr;
-    astar_pts[nxt_avail_pt].downc = downc;
+/*    astar_pts[nxt_avail_pt].downr = downr;
+    astar_pts[nxt_avail_pt].downc = downc; */
 
     nxt_avail_pt++;
 
@@ -288,8 +287,8 @@ int replace(			/* ele was in there */
     while (heap_run <= heap_size) {
 	now = heap_index[heap_run];
 	if (astar_pts[now].r == upr && astar_pts[now].c == upc) {
-	    astar_pts[now].downr = r;
-	    astar_pts[now].downc = c;
+	    /*astar_pts[now].downr = r;
+	    astar_pts[now].downc = c; */
 	    return 0;
 	}
 	heap_run++;

raster/r.watershed/ram/do_cum.c

@@ -6,8 +6,10 @@
 int do_cum(void)
 {
     SHORT r, c, dr, dc;
-    CELL is_swale, value, valued;
+    CELL is_swale, value, valued, aspect;
     int killer, threshold, count;
+    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
 
     G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
 
@@ -20,12 +22,18 @@ int do_cum(void)
 	G_percent(count++, do_points, 2);
 	killer = first_cum;
 	first_cum = astar_pts[killer].nxt;
-	if ((dr = astar_pts[killer].downr) > -1) {
-	    r = astar_pts[killer].r;
-	    c = astar_pts[killer].c;
-	    dc = astar_pts[killer].downc;
+	r = astar_pts[killer].r;
+	c = astar_pts[killer].c;
+	aspect = asp[SEG_INDEX(asp_seg, r, c)];
+	if (aspect) {
+	    dr = r + asp_r[ABS(aspect)];
+	    dc = c + asp_c[ABS(aspect)];
+	}
+	else
+	    dr = dc = -1;
+	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
 	    value = wat[SEG_INDEX(wat_seg, r, c)];
-	    if (ABS(value) >= threshold)
+	    if ((int)(ABS(value) + 0.5) >= threshold)
 		FLAG_SET(swale, r, c);
 	    valued = wat[SEG_INDEX(wat_seg, dr, dc)];
 	    if (value > 0) {
@@ -43,6 +51,13 @@ int do_cum(void)
 	    wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
 	    valued = ABS(valued) + 0.5;
 	    is_swale = FLAG_GET(swale, r, c);
+	    /* update asp for depression */
+	    if (is_swale && pit_flag) {
+		if (aspect > 0 && asp[SEG_INDEX(asp_seg, dr, dc)] == 0) {
+		    aspect = -aspect;
+		    asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+		}
+	    }
 	    if (is_swale || ((int)valued) >= threshold) {
 		FLAG_SET(swale, dr, dc);
 	    }
@@ -86,17 +101,19 @@ int do_cum_mfd(void)
     CELL is_swale;
     DCELL value, valued;
     int killer, threshold, count;
+    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
 
     /* MFD */
-    int mfd_cells, stream_cells, swale_cells, astar_not_set;
+    int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
     double *dist_to_nbr, *weight, sum_weight, max_weight;
-    int r_nbr, c_nbr, r_max, c_max, r_swale, c_swale, ct_dir, np_side;
+    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
     double dx, dy;
     CELL ele, ele_nbr, aspect, is_worked;
-    double prop, max_acc, max_swale;
+    double prop, max_acc;
     int workedon;
 
-    G_message(_("SECTION 3: Accumulating Surface Flow with MFD"));
+    G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
     G_debug(1, "MFD convergence factor set to %d.", c_fac);
 
     /* distances to neighbours */
@@ -129,22 +146,32 @@ int do_cum_mfd(void)
 	G_percent(count++, do_points, 2);
 	killer = first_cum;
 	first_cum = astar_pts[killer].nxt;
-	if ((dr = astar_pts[killer].downr) > -1) {
-	    r = astar_pts[killer].r;
-	    c = astar_pts[killer].c;
-	    dc = astar_pts[killer].downc;
+	r = astar_pts[killer].r;
+	c = astar_pts[killer].c;
+	aspect = asp[SEG_INDEX(asp_seg, r, c)];
+	if (aspect) {
+	    dr = r + asp_r[ABS(aspect)];
+	    dc = c + asp_c[ABS(aspect)];
+	}
+	else
+	    dr = dc = -1;
+	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = astar_pts[killer].downr) > -1) { */
 	    value = wat[SEG_INDEX(wat_seg, r, c)];
 	    valued = wat[SEG_INDEX(wat_seg, dr, dc)];
 
-	    /* start MFD */
+	    r_max = dr;
+	    c_max = dc;
 
 	    /* get weights */
 	    max_weight = 0;
 	    sum_weight = 0;
 	    np_side = -1;
 	    mfd_cells = 0;
+	    stream_cells = 0;
+	    swale_cells = 0;
 	    astar_not_set = 1;
 	    ele = alt[SEG_INDEX(alt_seg, r, c)];
+	    is_null = 0;
 	    /* this loop is needed to get the sum of weights */
 	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 		/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -154,19 +181,37 @@ int do_cum_mfd(void)
 		/* check that neighbour is within region */
 		if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 		    c_nbr < ncols) {
+
+		    /* check for swale or stream cells */
+		    is_swale = FLAG_GET(swale, r_nbr, c_nbr);
+		    if (is_swale)
+			swale_cells++;
+		    valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
+		    if ((ABS(valued) + 0.5) >= threshold)
+			stream_cells++;
+
 		    is_worked = FLAG_GET(worked, r_nbr, c_nbr);
 		    if (is_worked == 0) {
 			ele_nbr = alt[SEG_INDEX(alt_seg, r_nbr, c_nbr)];
-			if (ele_nbr < ele) {
-			    weight[ct_dir] =
-				mfd_pow(((ele -
-					  ele_nbr) / dist_to_nbr[ct_dir]),
-					c_fac);
+			is_null = G_is_c_null_value(&ele_nbr);
+			if (!is_null && ele_nbr <= ele) {
+			    if (ele_nbr < ele) {
+				weight[ct_dir] =
+				    mfd_pow(((ele -
+					      ele_nbr) / dist_to_nbr[ct_dir]),
+					    c_fac);
+			    }
+			    if (ele_nbr == ele) {
+				weight[ct_dir] =
+				    mfd_pow((0.5 / dist_to_nbr[ct_dir]),
+					    c_fac);
+			    }
 			    sum_weight += weight[ct_dir];
 			    mfd_cells++;
 
-			    if (weight[ct_dir] > max_weight)
+			    if (weight[ct_dir] > max_weight) {
 				max_weight = weight[ct_dir];
+			    }
 
 			    if (dr == r_nbr && dc == c_nbr) {
 				astar_not_set = 0;
@@ -181,7 +226,6 @@ int do_cum_mfd(void)
 	    /* honour A * path 
 	     * mfd_cells == 0: fine, SFD along A * path
 	     * mfd_cells == 1 && astar_not_set == 0: fine, SFD along A * path
-	     * mfd_cells > 1 && astar_not_set == 0: MFD, set A * path to max_weight
 	     * mfd_cells > 0 && astar_not_set == 1: A * path not included, add to mfd_cells
 	     */
 
@@ -192,21 +236,8 @@ int do_cum_mfd(void)
 		weight[np_side] = max_weight;
 	    }
 
-	    /* MFD, A * path included, set A * path to max_weight */
-	    if (mfd_cells > 1 && astar_not_set == 0) {
-		sum_weight = sum_weight - weight[np_side] + max_weight;
-		weight[np_side] = max_weight;
-	    }
-
 	    /* set flow accumulation for neighbours */
-	    stream_cells = 0;
-	    swale_cells = 0;
 	    max_acc = -1;
-	    max_swale = -1;
-	    r_max = dr;
-	    c_max = dc;
-	    r_swale = dr;
-	    c_swale = dc;
 
 	    if (mfd_cells > 1) {
 		prop = 0.0;
@@ -219,11 +250,11 @@ int do_cum_mfd(void)
 		    if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 			c_nbr < ncols && weight[ct_dir] > -0.5) {
 			is_worked = FLAG_GET(worked, r_nbr, c_nbr);
-			is_swale = FLAG_GET(swale, r_nbr, c_nbr);
 			if (is_worked == 0) {
 
 			    weight[ct_dir] = weight[ct_dir] / sum_weight;
-			    prop += weight[ct_dir];	/* check everything sums up to 1.0 */
+			    /* check everything sums up to 1.0 */
+			    prop += weight[ct_dir];
 
 			    valued = wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)];
 			    if (value > 0) {
@@ -240,39 +271,27 @@ int do_cum_mfd(void)
 			    }
 			    wat[SEG_INDEX(wat_seg, r_nbr, c_nbr)] = valued;
 
-			    if ((ABS(valued) + 0.5) >= threshold)
-				stream_cells++;
-			    if (ABS(valued) > max_acc) {
+			    /* get main drainage direction */
+			    if (ABS(valued) >= max_acc) {
 				max_acc = ABS(valued);
 				r_max = r_nbr;
 				c_max = c_nbr;
 			    }
-			    /* assist in stream merging */
-			    if (is_swale && ABS(valued) > max_swale) {
-				max_swale = ABS(valued);
-				r_swale = r_nbr;
-				c_swale = c_nbr;
-			    }
 			}
 			else if (ct_dir == np_side) {
-			    workedon++;	/* check for consistency with A * path */
+			    /* check for consistency with A * path */
+			    workedon++;
 			}
-			if (is_swale)
-			    swale_cells++;
 		    }
 		}
 		if (ABS(prop - 1.0) > 5E-6f) {
 		    G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
 			      prop);
 		}
-		valued = wat[SEG_INDEX(wat_seg, dr, dc)];
-		if (max_swale > -0.5) {
-		    r_max = r_swale;
-		    c_max = c_swale;
-		}
 	    }
 
 	    if (mfd_cells < 2) {
+		valued = wat[SEG_INDEX(wat_seg, dr, dc)];
 		if (value > 0) {
 		    if (valued > 0)
 			valued += value;
@@ -288,14 +307,27 @@ int do_cum_mfd(void)
 		wat[SEG_INDEX(wat_seg, dr, dc)] = valued;
 	    }
 
-	    /* MFD finished */
-
-	    valued = ABS(valued) + 0.5;
+	    /* update asp */
+	    if (dr != r_max || dc != c_max) {
+		aspect = drain[r - r_max + 1][c - c_max + 1];
+		if (asp[SEG_INDEX(asp_seg, r, c)] < 0)
+		    aspect = -aspect;
+		asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+	    }
 	    is_swale = FLAG_GET(swale, r, c);
-	    /* start new stream: only works with A * path */
-	    if (!is_swale && (int)valued >= threshold && stream_cells < 2 &&
+	    /* update asp for depression */
+	    if (is_swale && pit_flag) {
+		if (aspect > 0 && asp[SEG_INDEX(asp_seg, r_max, c_max)] == 0) {
+		    aspect = -aspect;
+		    asp[SEG_INDEX(asp_seg, r, c)] = aspect;
+		}
+	    }
+	    /* start new stream */
+	    value = ABS(value) + 0.5;
+	    if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
 		swale_cells < 1) {
-		FLAG_SET(swale, dr, dc);
+		FLAG_SET(swale, r, c);
+		is_swale = 1;
 	    }
 	    /* continue stream */
 	    if (is_swale) {
@@ -305,26 +337,21 @@ int do_cum_mfd(void)
 		if (er_flag && !is_swale)
 		    slope_length(r, c, r_max, c_max);
 	    }
-	    /* update asp */
-	    if (dr != r_max || dc != c_max) {
-		aspect = drain[r - r_max + 1][c - c_max + 1];
-		if (asp[SEG_INDEX(asp_seg, r, c)] < 0)
-		    aspect = -aspect;
-		asp[SEG_INDEX(asp_seg, r, c)] = aspect;
-
-	    }
 	    FLAG_SET(worked, r, c);
 	}
     }
     G_percent(count, do_points, 1);	/* finish it */
     if (workedon)
-	G_warning(_("MFD: A * path already processed when distributing flow: %d of %d"),
+	G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
 		  workedon, do_points);
 
     G_free(astar_pts);
 
     flag_destroy(worked);
 
+    G_free(dist_to_nbr);
+    G_free(weight);
+
     return 0;
 }
 

raster/r.watershed/ram/init_vars.c

@@ -8,7 +8,7 @@
 int init_vars(int argc, char *argv[])
 {
     SHORT r, c;
-    CELL *buf, alt_value, wat_value, asp_value;
+    CELL *buf, alt_value, wat_value, asp_value, block_value;
     int fd, index;
     char MASK_flag;
 
@@ -121,8 +121,10 @@ int init_vars(int argc, char *argv[])
     buf = G_allocate_cell_buf();
     alt =
 	(CELL *) G_malloc(sizeof(CELL) * size_array(&alt_seg, nrows, ncols));
-    r_h =
-	(CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
+    if (er_flag) {
+	r_h =
+	    (CELL *) G_malloc(sizeof(CELL) * size_array(&r_h_seg, nrows, ncols));
+    }
 
     fd = G_open_cell_old(ele_name, "");
     if (fd < 0) {
@@ -133,18 +135,31 @@ int init_vars(int argc, char *argv[])
     in_list = flag_create(nrows, ncols);
     worked = flag_create(nrows, ncols);
 
+    MASK_flag = 0;
+    do_points = nrows * ncols;
     for (r = 0; r < nrows; r++) {
 	G_get_c_raster_row(fd, buf, r);
 	for (c = 0; c < ncols; c++) {
 	    index = SEG_INDEX(alt_seg, r, c);
-	    alt[index] = r_h[index] = buf[c];
+	    alt_value = alt[index] = buf[c];
+	    if (er_flag) {
+		r_h[index] = buf[c];
+	    }
 	    /* all flags need to be manually set to zero */
-	    flag_unset(swale, r, c);	
-	    flag_unset(in_list, r, c);	
-	    flag_unset(worked, r, c);	
+	    flag_unset(swale, r, c);
+	    flag_unset(in_list, r, c);
+	    flag_unset(worked, r, c);
+	    /* check for masked and NULL cells here */
+	    if (G_is_c_null_value(&alt_value)) {
+		FLAG_SET(worked, r, c);
+		FLAG_SET(in_list, r, c);
+		do_points--;
+	    }
 	}
     }
     G_close_cell(fd);
+    if (do_points < nrows * ncols)
+	MASK_flag = 1;
     wat =
 	(DCELL *) G_malloc(sizeof(DCELL) *
 			   size_array(&wat_seg, nrows, ncols));
@@ -157,20 +172,36 @@ int init_vars(int argc, char *argv[])
 	for (r = 0; r < nrows; r++) {
 	    G_get_c_raster_row(fd, buf, r);
 	    for (c = 0; c < ncols; c++) {
-		wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
+		if (MASK_flag) {
+		    index = FLAG_GET(worked, r, c);
+		    if (!index)
+			wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
+		    else
+			wat[SEG_INDEX(wat_seg, r, c)] = 0.0;
+		}
+		else
+		    wat[SEG_INDEX(wat_seg, r, c)] = buf[c];
 	    }
 	}
 	G_close_cell(fd);
     }
     else {
 	for (r = 0; r < nrows; r++) {
-	    for (c = 0; c < ncols; c++)
-		wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+	    for (c = 0; c < ncols; c++) {
+		if (MASK_flag) {
+		    index = FLAG_GET(worked, r, c);
+		    if (!index)
+			wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+		}
+		else
+		    wat[SEG_INDEX(wat_seg, r, c)] = 1.0;
+	    }
 	}
     }
     asp =
 	(CELL *) G_malloc(size_array(&asp_seg, nrows, ncols) * sizeof(CELL));
 
+    /* depression: drainage direction will be set to zero later */
     if (pit_flag) {
 	fd = G_open_cell_old(pit_name, "");
 	if (fd < 0) {
@@ -179,12 +210,15 @@ int init_vars(int argc, char *argv[])
 	for (r = 0; r < nrows; r++) {
 	    G_get_c_raster_row(fd, buf, r);
 	    for (c = 0; c < ncols; c++) {
-		asp[SEG_INDEX(asp_seg, r, c)] = buf[c];
+		asp_value = buf[c];
+		if (!G_is_c_null_value(&asp_value))
+		    asp[SEG_INDEX(asp_seg, r, c)] = 1;
 	    }
 	}
 	G_close_cell(fd);
     }
 
+    /* this is also creating streams... */
     if (ob_flag) {
 	fd = G_open_cell_old(ob_name, "");
 	if (fd < 0) {
@@ -193,7 +227,8 @@ int init_vars(int argc, char *argv[])
 	for (r = 0; r < nrows; r++) {
 	    G_get_c_raster_row(fd, buf, r);
 	    for (c = 0; c < ncols; c++) {
-		if (buf[c])
+		block_value = buf[c];
+		if (!G_is_c_null_value(&block_value))
 		    FLAG_SET(swale, r, c);
 	    }
 	}
@@ -206,32 +241,12 @@ int init_vars(int argc, char *argv[])
 	}
     }
 
-    MASK_flag = 0;
-    do_points = nrows * ncols;
-    if (NULL != G_find_file("cell", "MASK", G_mapset())) {
-	MASK_flag = 1;
-	if ((fd = G_open_cell_old("MASK", G_mapset())) < 0) {
-	    G_fatal_error(_("unable to open MASK"));
-	}
-	else {
-	    for (r = 0; r < nrows; r++) {
-		G_get_c_raster_row_nomask(fd, buf, r);
-		for (c = 0; c < ncols; c++) {
-		    if (!buf[c]) {
-			FLAG_SET(worked, r, c);
-			FLAG_SET(in_list, r, c);
-			do_points--;
-		    }
-		}
-	    }
-	    G_close_cell(fd);
-	}
+    /* RUSLE: LS and/or S factor */
+    if (er_flag) {
+	s_l =
+	    (double *)G_malloc(size_array(&s_l_seg, nrows, ncols) *
+			       sizeof(double));
     }
-    s_l =
-	(double *)G_malloc(size_array(&s_l_seg, nrows, ncols) *
-			   sizeof(double));
-    /* astar_pts = (POINT *) G_malloc (nrows * ncols * sizeof (POINT)); */
-    astar_pts = (POINT *) G_malloc(do_points * sizeof(POINT));
 
     if (sg_flag) {
 	s_g =
@@ -244,6 +259,8 @@ int init_vars(int argc, char *argv[])
 			       sizeof(double));
     }
 
+    astar_pts = (POINT *) G_malloc(do_points * sizeof(POINT));
+
     /* heap_index will track astar_pts in ternary min-heap */
     /* heap_index is one-based */
     heap_index = (int *)G_malloc((do_points + 1) * sizeof(int));
@@ -262,14 +279,20 @@ int init_vars(int argc, char *argv[])
 		    wat[SEG_INDEX(wat_seg, r, c)] = 0;
 		}
 		else {
-		    s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+		    if (er_flag)
+			s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
 		    asp_value = asp[SEG_INDEX(asp_seg, r, c)];
 		    if (r == 0 || c == 0 || r == nrows - 1 ||
 			c == ncols - 1 || asp_value != 0) {
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
 			    wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
-			if (r == 0)
+			/* set depression */
+			if (asp_value) {
+			    asp_value = 0;
+			    wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+			}
+			else if (r == 0)
 			    asp_value = -2;
 			else if (c == 0)
 			    asp_value = -4;
@@ -277,15 +300,13 @@ int init_vars(int argc, char *argv[])
 			    asp_value = -6;
 			else if (c == ncols - 1)
 			    asp_value = -8;
-			else
-			    asp_value = -1;
 			asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 		    }
 		    else if (FLAG_GET(worked, r - 1, c)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -2;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -293,7 +314,7 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else if (FLAG_GET(worked, r + 1, c)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -6;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -301,7 +322,7 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else if (FLAG_GET(worked, r, c - 1)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -4;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -309,7 +330,7 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else if (FLAG_GET(worked, r, c + 1)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -8;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -317,7 +338,7 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r - 1, c - 1)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -3;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -325,7 +346,7 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r - 1, c + 1)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -1;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -333,7 +354,7 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r + 1, c - 1)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -5;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -341,7 +362,7 @@ int init_vars(int argc, char *argv[])
 		    }
 		    else if (sides == 8 && FLAG_GET(worked, r + 1, c + 1)) {
 			alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp[SEG_INDEX(asp_seg, r, c)] = -7;
 			wat_value = wat[SEG_INDEX(wat_seg, r, c)];
 			if (wat_value > 0)
@@ -355,7 +376,8 @@ int init_vars(int argc, char *argv[])
 	for (r = 0; r < nrows; r++) {
 	    G_percent(r, nrows, 3);
 	    for (c = 0; c < ncols; c++) {
-		s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
+		if (er_flag)
+		    s_l[SEG_INDEX(s_l_seg, r, c)] = half_res;
 		asp_value = asp[SEG_INDEX(asp_seg, r, c)];
 		if (r == 0 || c == 0 || r == nrows - 1 ||
 		    c == ncols - 1 || asp_value != 0) {
@@ -363,7 +385,12 @@ int init_vars(int argc, char *argv[])
 		    if (wat_value > 0) {
 			wat[SEG_INDEX(wat_seg, r, c)] = -wat_value;
 		    }
-		    if (r == 0)
+		    /* set depression */
+		    if (asp_value) {
+			asp_value = 0;
+			wat[SEG_INDEX(wat_seg, r, c)] = ABS(wat_value);
+		    }
+		    else if (r == 0)
 			asp_value = -2;
 		    else if (c == 0)
 			asp_value = -4;
@@ -371,17 +398,15 @@ int init_vars(int argc, char *argv[])
 			asp_value = -6;
 		    else if (c == ncols - 1)
 			asp_value = -8;
-		    else
-			asp_value = -1;
 		    asp[SEG_INDEX(asp_seg, r, c)] = asp_value;
 		    alt_value = alt[SEG_INDEX(alt_seg, r, c)];
-		    add_pt(r, c, -1, -1, alt_value, alt_value);
+		    add_pt(r, c, alt_value, alt_value);
 		}
 	    }
 	}
     }
 
-    G_percent(r, nrows, 3);	/* finish it */
+    G_percent(r, nrows, 1);	/* finish it */
 
     return 0;
 }

raster/r.watershed/ram/sg_factor.c

@@ -9,7 +9,7 @@ int sg_factor(void)
     CELL low_elev, hih_elev;
     double height, length, S, sin_theta;
 
-    G_message(_("SECTION 4: Length Slope determination."));
+    G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
 
     for (r = 0; r < nrows; r++) {
 	G_percent(r, nrows, 3);
@@ -38,7 +38,7 @@ int sg_factor(void)
 	    }
 	}
     }
-    G_percent(r, nrows, 3);	/* finish it */
+    G_percent(nrows, nrows, 1);	/* finish it */
 
     if (ril_flag) {
 	G_free(ril_buf);

raster/r.watershed/ram/slope_len.c

@@ -20,7 +20,7 @@ int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
 	    if (asp_value == 2 || asp_value == 6)
 		res = window.ns_res;
 	    else		/* asp_value == 4, 8, -2, -4, -6, or -8 */
-		res = diag;
+		res = diag;     /* how can res be diag with sides == 4??? */
 	}
 	else {			/* c == dc */
 

raster/r.watershed/seg/Gwater.h

@@ -31,7 +31,7 @@
 
 #define POINT       struct points
 POINT {
-    SHORT r, c, downr, downc;
+    SHORT r, c; /* , downr, downc */
     int nxt;
 };
 
@@ -87,7 +87,7 @@ CELL def_basin(int, int, CELL, double, CELL);
 
 /* do_astar.c */
 int do_astar(void);
-int add_pt(SHORT, SHORT, SHORT, SHORT, CELL, CELL);
+int add_pt(SHORT, SHORT, CELL, CELL);
 int drop_pt(void);
 int sift_up(int, CELL);
 double get_slope(SHORT, SHORT, SHORT, SHORT, CELL, CELL);

raster/r.watershed/seg/close_maps.c

@@ -1,69 +1,110 @@
 #include "Gwater.h"
 #include <unistd.h>
+#include <grass/glocale.h>
 
 int close_maps(void)
 {
-    struct Colors colors, logcolors;
+    struct Colors colors;
     int r, c;
     CELL is_swale;
-    double dvalue;
+    DCELL *dbuf = NULL;
+    int fd;
     struct FPRange accRange;
     DCELL min, max;
     DCELL clr_min, clr_max;
+    DCELL sum, sum_sqr, stddev, lstddev, dvalue;
 
     dseg_close(&slp);
     cseg_close(&alt);
     if (wat_flag) {
 	dseg_write_cellfile(&wat, wat_name);
 
+	/* get standard deviation */
+	fd = G_open_cell_old(wat_name, "");
+	if (fd < 0) {
+	    G_fatal_error(_("unable to open flow accumulation map layer"));
+	}
+
+	sum = sum_sqr = stddev = 0.0;
+	dbuf = G_allocate_d_raster_buf();
+	for (r = 0; r < nrows; r++) {
+	    G_get_d_raster_row(fd, dbuf, r);
+	    for (c = 0; c < ncols; c++) {
+		dvalue = dbuf[c];
+		if (G_is_d_null_value(&dvalue) == 0 && dvalue) {
+		    dvalue = ABS(dvalue);
+		    sum += dvalue;
+		    sum_sqr += dvalue * dvalue;
+		}
+	    }
+	}
+
+	stddev = sqrt((sum_sqr - (sum + sum / do_points)) / (do_points - 1));
+	G_debug(1, "stddev: %f", stddev);
+
 	/* set nice color rules: yellow, green, cyan, blue, black */
+	/* start with white to get more detail? NULL cells are white by default, may be confusing */
+
+	lstddev = log(stddev);
+
 	G_read_fp_range(wat_name, this_mapset, &accRange);
 	min = max = 0;
 	G_get_fp_range_min_max(&accRange, &min, &max);
 
 	G_init_colors(&colors);
-	if (ABS(min) > max) {
-	    clr_min = 1;
-	    clr_max = 0.25 * max;
-	    G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0, 255,
-				      0, &colors);
-	    clr_min = 0.25 * max;
-	    clr_max = 0.5 * max;
-	    G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0, 255,
-				      255, &colors);
-	    clr_min = 0.5 * max;
-	    clr_max = 0.75 * max;
-	    G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0, 0,
-				      255, &colors);
-	    clr_min = 0.75 * max;
-	    clr_max = max;
-	    G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0, 0,
-				      &colors);
-	    max = -max;
+
+	if (min < 0) {
+	    if (min < (-stddev - 1)) {
+		clr_min = min;
+		clr_max = -stddev - 1;
+		G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+					  0, 0, &colors);
+	    }
+	    clr_min = -stddev - 1.;
+	    clr_max = -1. * exp(lstddev * 0.75);
+	    G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+				      0, 255, &colors);
+	    clr_min = clr_max;
+	    clr_max = -1. * exp(lstddev * 0.5);
+	    G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0,
+				      255, 255, &colors);
+	    clr_min = clr_max;
+	    clr_max = -1. * exp(lstddev * 0.35);
+	    G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+				      255, 0, &colors);
+	    clr_min = clr_max;
+	    clr_max = -1.;
+	    G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 255,
+				      255, 0, &colors);
 	}
-	else {
-	    min = ABS(min);
-	    clr_min = 1;
-	    clr_max = 0.25 * min;
-	    G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0, 255,
-				      0, &colors);
-	    clr_min = 0.25 * min;
-	    clr_max = 0.5 * min;
-	    G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0, 255,
-				      255, &colors);
-	    clr_min = 0.5 * min;
-	    clr_max = 0.75 * min;
-	    G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0, 0,
-				      255, &colors);
-	    clr_min = 0.75 * min;
-	    clr_max = min;
-	    G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0, 0,
+	clr_min = -1.;
+	clr_max = 1.;
+	G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 255,
+				  255, 0, &colors);
+	clr_min = 1;
+	clr_max = exp(lstddev * 0.35);
+	G_add_d_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
+				  255, 0, &colors);
+	clr_min = clr_max;
+	clr_max = exp(lstddev * 0.5);
+	G_add_d_raster_color_rule(&clr_min, 0, 255, 0, &clr_max, 0,
+				  255, 255, &colors);
+	clr_min = clr_max;
+	clr_max = exp(lstddev * 0.75);
+	G_add_d_raster_color_rule(&clr_min, 0, 255, 255, &clr_max, 0,
+				  0, 255, &colors);
+	clr_min = clr_max;
+	clr_max = stddev + 1.;
+	G_add_d_raster_color_rule(&clr_min, 0, 0, 255, &clr_max, 0, 0,
+				  0, &colors);
+
+	if (max > 0 && max > stddev + 1) {
+	    clr_min = stddev + 1;
+	    clr_max = max;
+	    G_add_d_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0, 0, 0,
 				      &colors);
 	}
-	G_abs_log_colors(&logcolors, &colors, 5);
-	G_add_d_raster_color_rule(&min, 0, 0, 0, &max, 0, 0, 0, &logcolors);
-	G_write_colors(wat_name, this_mapset, &logcolors);
-
+	G_write_colors(wat_name, this_mapset, &colors);
     }
     if (asp_flag) {
 	cseg_write_cellfile(&asp, asp_name);

raster/r.watershed/seg/do_astar.c

@@ -46,7 +46,8 @@ int do_astar(void)
 	c = point.c;
 
 	bseg_put(&worked, &one, r, c);
-	cseg_get(&alt, &alt_val, r, c);
+	/* cseg_get(&alt, &alt_val, r, c); */
+	alt_val = heap_pos.ele;
 
 	/* check all neighbours, breadth first search */
 	for (ct_dir = 0; ct_dir < sides; ct_dir++) {
@@ -60,26 +61,27 @@ int do_astar(void)
 		bseg_get(&in_list, &in_val, upr, upc);
 		if (in_val == 0) {
 		    cseg_get(&alt, &alt_up, upr, upc);
-		    add_pt(upr, upc, r, c, alt_up, alt_val);
+		    add_pt(upr, upc, alt_up, alt_val);
 		    /* flow direction is set here */
 		    drain_val = drain[upr - r + 1][upc - c + 1];
 		    cseg_put(&asp, &drain_val, upr, upc);
 		}
 		else {
 		    /* check if neighbour has not been worked on,
-		     * update values for asp, wat and slp */
+		     * update values for asp and wat */
 		    bseg_get(&worked, &work_val, upr, upc);
 		    if (!work_val) {
 			cseg_get(&asp, &asp_up, upr, upc);
-			if (asp_up < -1) {
+			if (asp_up < 0) {
 			    drain_val = drain[upr - r + 1][upc - c + 1];
 			    cseg_put(&asp, &drain_val, upr, upc);
 			    dseg_get(&wat, &wat_val, r, c);
-			    if (wat_val > 0)
+			    if (wat_val > 0) {
 				wat_val = -wat_val;
-			    dseg_put(&wat, &wat_val, r, c);
-			    cseg_get(&alt, &alt_up, upr, upc);
-			    replace(upr, upc, r, c);	/* alt_up used to be */
+				dseg_put(&wat, &wat_val, r, c);
+			    }
+			    /* cseg_get(&alt, &alt_up, upr, upc); */
+			    /* replace(upr, upc, r, c); */	/* alt_up used to be */
 			    /* slope = get_slope (upr, upc, r, c, alt_up, alt_val);
 			       dseg_put (&slp, &slope, upr, upc); */
 			}
@@ -100,7 +102,7 @@ int do_astar(void)
 }
 
 /* new add point routine for min heap */
-int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
+int add_pt(SHORT r, SHORT c, CELL ele, CELL downe)
 {
     POINT point;
     HEAP heap_pos;
@@ -124,11 +126,13 @@ int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 
     point.r = r;
     point.c = c;
-    point.downr = downr;
-    point.downc = downc;
+/*    point.downr = downr;
+    point.downc = downc; */
 
     seg_put(&astar_pts, (char *)&point, 0, nxt_avail_pt);
 
+    /* cseg_put(&pnt_index, &nxt_avail_pt, r, c); */
+
     nxt_avail_pt++;
 
     /* sift up: move new point towards top of heap */
@@ -142,8 +146,8 @@ int add_pt(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
 int drop_pt(void)
 {
     int child, childr, parent;
-    int childp, childrp;
-    CELL ele, eler;
+    int childp; /* , childrp */
+    CELL ele; /* , eler */
     int i;
     HEAP heap_pos;
 
@@ -170,23 +174,23 @@ int drop_pt(void)
 	ele = heap_pos.ele;
 	if (child < heap_size) {
 	    childr = child + 1;
-	    i = 1;
-	    while (childr <= heap_size && i < 3) {
+	    i = child + 3;
+	    while (childr <= heap_size && childr < i) {
 		seg_get(&heap_index, (char *)&heap_pos, 0, childr);
-		childrp = heap_pos.point;
-		eler = heap_pos.ele;
-		if (eler < ele) {
+		/* childrp = heap_pos.point;
+		eler = heap_pos.ele; */
+		if (heap_pos.ele < ele) {
 		    child = childr;
-		    childp = childrp;
-		    ele = eler;
+		    childp = heap_pos.point;
+		    ele = heap_pos.ele;
 		}
 		/* make sure we get the oldest child */
-		else if (ele == eler && childp > childrp) {
+		else if (ele == heap_pos.ele && childp > heap_pos.point) {
 		    child = childr;
-		    childp = childrp;
+		    childp = heap_pos.point;
 		}
 		childr++;
-		i++;
+		/* i++; */
 	    }
 	}
 
@@ -220,8 +224,8 @@ int drop_pt(void)
 /* standard sift-up routine for d-ary min heap */
 int sift_up(int start, CELL ele)
 {
-    int parent, parentp, child, childp;
-    CELL elep;
+    int parent, child, childp; /* parentp, */
+    /* CELL elep; */
     HEAP heap_pos;
 
     child = start;
@@ -231,11 +235,11 @@ int sift_up(int start, CELL ele)
     while (child > 1) {
 	parent = GET_PARENT(child);
 	seg_get(&heap_index, (char *)&heap_pos, 0, parent);
-	parentp = heap_pos.point;
-	elep = heap_pos.ele;
+	/* parentp = heap_pos.point;
+	elep = heap_pos.ele; */
 
 	/* parent ele higher */
-	if (elep > ele) {
+	if (heap_pos.ele > ele) {
 
 	    /* push parent point down */
 	    seg_put(&heap_index, (char *)&heap_pos, 0, child);
@@ -243,7 +247,7 @@ int sift_up(int start, CELL ele)
 
 	}
 	/* same ele, but parent is younger */
-	else if (elep == ele && parentp > childp) {
+	else if (heap_pos.ele == ele && heap_pos.point > childp) {
 
 	    /* push parent point down */
 	    seg_put(&heap_index, (char *)&heap_pos, 0, child);
@@ -280,28 +284,23 @@ get_slope(SHORT r, SHORT c, SHORT downr, SHORT downc, CELL ele, CELL downe)
     return (slope);
 }
 
-int replace(			/* ele was in there */
-	       SHORT upr, SHORT upc, SHORT r, SHORT c)
-/* CELL ele;  */
-{
-    int now, heap_run;
+int replace(SHORT upr, SHORT upc, SHORT r, SHORT c)
+{				/* ele was in there */
+    /* CELL ele;  */
+    int now;
     POINT point;
-    HEAP heap_pos;
 
-    heap_run = 0;
-
-    /* this is dumb, works for ram, for seg make new index pt_index[row][col] */
-    while (heap_run <= heap_size) {
-	seg_get(&heap_index, (char *)&heap_pos, 0, heap_run);
-	now = heap_pos.point;
-	seg_get(&astar_pts, (char *)&point, 0, now);
-	if (point.r == upr && point.c == upc) {
-	    point.downr = r;
-	    point.downc = c;
-	    seg_put(&astar_pts, (char *)&point, 0, now);
-	    return 0;
-	}
-	heap_run++;;
+    now = 0;
+
+    /* cseg_get(&pnt_index, &now, upr, upc); */
+    seg_get(&astar_pts, (char *)&point, 0, now);
+    if (point.r != upr || point.c != upc) {
+	G_warning("pnt_index incorrect!");
+	return 1;
     }
+/*    point.downr = r;
+    point.downc = c; */
+    seg_put(&astar_pts, (char *)&point, 0, now);
+
     return 0;
 }

raster/r.watershed/seg/do_cum.c

@@ -7,13 +7,14 @@
 int do_cum(void)
 {
     SHORT r, c, dr, dc;
-    CELL is_swale;
+    CELL is_swale, asp_val, asp_val_down;
     DCELL value, valued;
     POINT point;
     int killer, threshold, count;
+    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
 
     G_message(_("SECTION 3: Accumulating Surface Flow with SFD."));
-    G_debug(1, "MFD convergence factor set to %d.", c_fac);
 
     count = 0;
     if (bas_thres <= 0)
@@ -25,10 +26,16 @@ int do_cum(void)
 	killer = first_cum;
 	seg_get(&astar_pts, (char *)&point, 0, killer);
 	first_cum = point.nxt;
-	if ((dr = point.downr) > -1) {
-	    r = point.r;
-	    c = point.c;
-	    dc = point.downc;
+	r = point.r;
+	c = point.c;
+	cseg_get(&asp, &asp_val, r, c);
+	if (asp_val) {
+	    dr = r + asp_r[ABS(asp_val)];
+	    dc = c + asp_c[ABS(asp_val)];
+	}
+	else
+	    dr = dc = -1;
+	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = point.downr) > -1) { */
 	    dseg_get(&wat, &value, r, c);
 	    if (ABS(value) >= threshold)
 		bseg_put(&swale, &one, r, c);
@@ -47,6 +54,14 @@ int do_cum(void)
 	    }
 	    dseg_put(&wat, &valued, dr, dc);
 	    bseg_get(&swale, &is_swale, r, c);
+	    /* update asp for depression */
+	    if (is_swale && pit_flag) {
+		cseg_get(&asp, &asp_val_down, dr, dc);
+		if (asp_val > 0 && asp_val_down == 0) {
+		    asp_val = -asp_val;
+		    cseg_put(&asp, &asp_val, r, c);
+		}
+	    }
 	    if (is_swale || ABS(valued) >= threshold) {
 		bseg_put(&swale, &one, dr, dc);
 	    }
@@ -88,20 +103,23 @@ int do_cum_mfd(void)
 {
     int r, c, dr, dc;
     CELL is_swale;
-    DCELL value, valued;
+    DCELL value, valued, *wat_nbr;
     POINT point;
     int killer, threshold, count;
 
     /* MFD */
-    int mfd_cells, stream_cells, swale_cells, astar_not_set;
+    int mfd_cells, stream_cells, swale_cells, astar_not_set, is_null;
     double *dist_to_nbr, *weight, sum_weight, max_weight;
-    int r_nbr, c_nbr, r_max, c_max, r_swale, c_swale, ct_dir, np_side;
+    int r_nbr, c_nbr, r_max, c_max, ct_dir, np_side;
     double dx, dy;
-    CELL ele, ele_nbr, aspect, asp_val, is_worked;
-    double prop, max_acc, max_swale;
+    CELL ele, ele_nbr, asp_val, asp_val2, cvalue, *worked_nbr;
+    double prop, max_acc;
     int workedon;
+    SHORT asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
+    SHORT asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
 
     G_message(_("SECTION 3: Accumulating Surface Flow with MFD."));
+    G_debug(1, "MFD convergence factor set to %d.", c_fac);
 
     /* distances to neighbours */
     dist_to_nbr = (double *)G_malloc(sides * sizeof(double));
@@ -127,6 +145,9 @@ int do_cum_mfd(void)
 	}
     }
 
+    worked_nbr = (CELL *)G_malloc(sides * sizeof(CELL));
+    wat_nbr = (DCELL *)G_malloc(sides * sizeof(DCELL));
+
     workedon = 0;
 
     count = 0;
@@ -140,26 +161,32 @@ int do_cum_mfd(void)
 	killer = first_cum;
 	seg_get(&astar_pts, (char *)&point, 0, killer);
 	first_cum = point.nxt;
-	if ((dr = point.downr) > -1) {
-	    r = point.r;
-	    c = point.c;
-	    dc = point.downc;
+	r = point.r;
+	c = point.c;
+	cseg_get(&asp, &asp_val, r, c);
+	if (asp_val) {
+	    dr = r + asp_r[ABS(asp_val)];
+	    dc = c + asp_c[ABS(asp_val)];
+	}
+	else
+	    dr = dc = -1;
+	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) { /* if ((dr = point.downr) > -1) { */
+	    /* dc = point.downc; */
 	    dseg_get(&wat, &value, r, c);
-	    dseg_get(&wat, &valued, dr, dc);
 
-	    /* disabled for MFD */
-	    /* if (ABS(value) > threshold)
-	       bseg_put(&swale, &one, r, c); */
-
-	    /* start MFD */
+	    r_max = dr;
+	    c_max = dc;
 
 	    /* get weights */
 	    max_weight = 0;
 	    sum_weight = 0;
 	    np_side = -1;
 	    mfd_cells = 0;
+	    stream_cells = 0;
+	    swale_cells = 0;
 	    astar_not_set = 1;
 	    cseg_get(&alt, &ele, r, c);
+	    is_null = 0;
 	    /* this loop is needed to get the sum of weights */
 	    for (ct_dir = 0; ct_dir < sides; ct_dir++) {
 		/* get r, c (r_nbr, c_nbr) for neighbours */
@@ -169,19 +196,39 @@ int do_cum_mfd(void)
 		/* check that neighbour is within region */
 		if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 		    c_nbr < ncols) {
-		    bseg_get(&worked, &is_worked, r_nbr, c_nbr);
-		    if (is_worked == 0) {
+
+		    /* check for swale or stream cells */
+		    bseg_get(&swale, &is_swale, r_nbr, c_nbr);
+		    if (is_swale)
+			swale_cells++;
+		    dseg_get(&wat, &valued, r_nbr, c_nbr);
+		    wat_nbr[ct_dir] = valued;
+		    if ((ABS(wat_nbr[ct_dir]) + 0.5) >= threshold)
+			stream_cells++;
+
+		    bseg_get(&worked, &cvalue, r_nbr, c_nbr);
+		    worked_nbr[ct_dir] = cvalue;
+		    if (worked_nbr[ct_dir] == 0) {
 			cseg_get(&alt, &ele_nbr, r_nbr, c_nbr);
-			if (ele_nbr < ele) {
-			    weight[ct_dir] =
-				mfd_pow(((ele -
-					  ele_nbr) / dist_to_nbr[ct_dir]),
-					c_fac);
+			is_null = G_is_c_null_value(&ele_nbr);
+			if (!is_null && ele_nbr <= ele) {
+			    if (ele_nbr < ele) {
+				weight[ct_dir] =
+				    mfd_pow(((ele -
+					      ele_nbr) / dist_to_nbr[ct_dir]),
+					    c_fac);
+			    }
+			    if (ele_nbr == ele) {
+				weight[ct_dir] =
+				    mfd_pow((0.5 / dist_to_nbr[ct_dir]),
+					    c_fac);
+			    }
 			    sum_weight += weight[ct_dir];
 			    mfd_cells++;
 
-			    if (weight[ct_dir] > max_weight)
+			    if (weight[ct_dir] > max_weight) {
 				max_weight = weight[ct_dir];
+			    }
 
 			    if (dr == r_nbr && dc == c_nbr) {
 				astar_not_set = 0;
@@ -196,7 +243,6 @@ int do_cum_mfd(void)
 	    /* honour A * path 
 	     * mfd_cells == 0: fine, SFD along A * path
 	     * mfd_cells == 1 && astar_not_set == 0: fine, SFD along A * path
-	     * mfd_cells > 1 && astar_not_set == 0: MFD, set A * path to max_weight
 	     * mfd_cells > 0 && astar_not_set == 1: A * path not included, add to mfd_cells
 	     */
 
@@ -207,21 +253,8 @@ int do_cum_mfd(void)
 		weight[np_side] = max_weight;
 	    }
 
-	    /* MFD, A * path included, set A * path to max_weight */
-	    if (mfd_cells > 1 && astar_not_set == 0) {
-		sum_weight = sum_weight - weight[np_side] + max_weight;
-		weight[np_side] = max_weight;
-	    }
-
 	    /* set flow accumulation for neighbours */
-	    stream_cells = 0;
-	    swale_cells = 0;
 	    max_acc = -1;
-	    max_swale = -1;
-	    r_max = dr;
-	    c_max = dc;
-	    r_swale = dr;
-	    c_swale = dc;
 
 	    if (mfd_cells > 1) {
 		prop = 0.0;
@@ -233,44 +266,38 @@ int do_cum_mfd(void)
 		    /* check that neighbour is within region */
 		    if (r_nbr >= 0 && r_nbr < nrows && c_nbr >= 0 &&
 			c_nbr < ncols && weight[ct_dir] > -0.5) {
-			bseg_get(&worked, &is_worked, r_nbr, c_nbr);
-			bseg_get(&swale, &is_swale, r_nbr, c_nbr);
-			if (is_worked == 0) {
+			/* bseg_get(&worked, &is_worked, r_nbr, c_nbr); */
+			if (worked_nbr[ct_dir] == 0) {
 
 			    weight[ct_dir] = weight[ct_dir] / sum_weight;
-			    prop += weight[ct_dir];	/* check everything sums up to 1.0 */
+			    /* check everything sums up to 1.0 */
+			    prop += weight[ct_dir];
 
-			    dseg_get(&wat, &valued, r_nbr, c_nbr);
 			    if (value > 0) {
-				if (valued > 0)
-				    valued += value * weight[ct_dir];
+				if (wat_nbr[ct_dir] > 0)
+				    wat_nbr[ct_dir] += value * weight[ct_dir];
 				else
-				    valued -= value * weight[ct_dir];
+				    wat_nbr[ct_dir] -= value * weight[ct_dir];
 			    }
 			    else {
-				if (valued < 0)
-				    valued += value * weight[ct_dir];
+				if (wat_nbr[ct_dir] < 0)
+				    wat_nbr[ct_dir] += value * weight[ct_dir];
 				else
-				    valued = value * weight[ct_dir] - valued;
+				    wat_nbr[ct_dir] = value * weight[ct_dir] - wat_nbr[ct_dir];
 			    }
+			    valued = wat_nbr[ct_dir];
 			    dseg_put(&wat, &valued, r_nbr, c_nbr);
 
-			    if ((ABS(valued) + 0.5) >= threshold)
-				stream_cells++;
-			    if (ABS(valued) > max_acc) {
-				max_acc = ABS(valued);
+			    /* get main drainage direction */
+			    if (ABS(wat_nbr[ct_dir]) >= max_acc) {
+				max_acc = ABS(wat_nbr[ct_dir]);
 				r_max = r_nbr;
 				c_max = c_nbr;
 			    }
-			    /* assist in stream merging */
-			    if (is_swale && ABS(valued) > max_swale) {
-				max_swale = ABS(valued);
-				r_swale = r_nbr;
-				c_swale = c_nbr;
-			    }
 			}
 			else if (ct_dir == np_side) {
-			    workedon++;	/* check for consistency with A * path */
+			    /* check for consistency with A * path */
+			    workedon++;
 			}
 		    }
 		}
@@ -278,14 +305,10 @@ int do_cum_mfd(void)
 		    G_warning(_("MFD: cumulative proportion of flow distribution not 1.0 but %f"),
 			      prop);
 		}
-		dseg_get(&wat, &valued, dr, dc);
-		if (max_swale > -0.5) {
-		    r_max = r_swale;
-		    c_max = c_swale;
-		}
 	    }
 
 	    if (mfd_cells < 2) {
+		dseg_get(&wat, &valued, dr, dc);
 		if (value > 0) {
 		    if (valued > 0)
 			valued += value;
@@ -301,14 +324,30 @@ int do_cum_mfd(void)
 		dseg_put(&wat, &valued, dr, dc);
 	    }
 
-	    /* MFD finished */
+	    /* update asp */
+	    if (dr != r_max || dc != c_max) {
+		asp_val2 = drain[r - r_max + 1][c - c_max + 1];
+		/* cseg_get(&asp, &asp_val, r, c); */
+		if (asp_val < 0)
+		    asp_val2 = -asp_val2;
+		cseg_put(&asp, &asp_val2, r, c);
 
-	    valued = ABS(valued) + 0.5;
+	    }
+	    /* update asp for depression */
 	    bseg_get(&swale, &is_swale, r, c);
-	    /* start new stream: only works with A * path */
-	    if (!is_swale && (int)valued >= threshold && stream_cells < 2 &&
+	    if (is_swale && pit_flag) {
+		cseg_get(&asp, &asp_val2, r_max, c_max);
+		if (asp_val > 0 && asp_val2 == 0) {
+		    asp_val = -asp_val;
+		    cseg_put(&asp, &asp_val, r, c);
+		}
+	    }
+	    /* start new stream */
+	    value = ABS(value) + 0.5;
+	    if (!is_swale && (int)value >= threshold && stream_cells < 1 &&
 		swale_cells < 1) {
-		bseg_put(&swale, &one, dr, dc);
+		bseg_put(&swale, &one, r, c);
+		is_swale = 1;
 	    }
 	    /* continue stream */
 	    if (is_swale) {
@@ -318,26 +357,22 @@ int do_cum_mfd(void)
 		if (er_flag && !is_swale)
 		    slope_length(r, c, r_max, c_max);
 	    }
-	    /* update asp */
-	    if (dr != r_max || dc != c_max) {
-		aspect = drain[r - r_max + 1][c - c_max + 1];
-		cseg_get(&asp, &asp_val, r, c);
-		if (asp_val < 0)
-		    aspect = -aspect;
-		cseg_put(&asp, &aspect, r, c);
-
-	    }
 	    bseg_put(&worked, &one, r, c);
 	}
     }
     G_percent(count, do_points, 1);	/* finish it */
     if (workedon)
-	G_warning(_("MFD: A * path already processed when distributing flow: %d of %d"),
+	G_warning(_("MFD: A * path already processed when distributing flow: %d of %d cells"),
 		  workedon, do_points);
 
     seg_close(&astar_pts);
 
     bseg_close(&worked);
+    
+    G_free(dist_to_nbr);
+    G_free(weight);
+    G_free(wat_nbr);
+    G_free(worked_nbr);
 
     return 0;
 }

raster/r.watershed/seg/dseg_get.c

@@ -4,7 +4,7 @@
 
 int dseg_get(DSEG * dseg, double *value, int row, int col)
 {
-    if (segment_get(&(dseg->seg), (CELL *) value, row, col) < 0) {
+    if (segment_get(&(dseg->seg), (DCELL *) value, row, col) < 0) {
 	G_warning("dseg_get(): could not read segment file");
 	return -1;
     }

raster/r.watershed/seg/dseg_put.c

@@ -4,7 +4,7 @@
 
 int dseg_put(DSEG * dseg, double *value, int row, int col)
 {
-    if (segment_put(&(dseg->seg), (CELL *) value, row, col) < 0) {
+    if (segment_put(&(dseg->seg), (DCELL *) value, row, col) < 0) {
 	G_warning("dseg_put(): could not write segment file");
 	return -1;
     }

raster/r.watershed/seg/dseg_read.c

@@ -6,10 +6,9 @@ static char *me = "dseg_read_cell";
 
 int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
 {
-    int row, col, nrows, ncols;
+    int row, nrows, ncols;
     int map_fd;
     char msg[100];
-    CELL *buffer;
     double *dbuffer;
 
     dseg->name = NULL;
@@ -24,11 +23,9 @@ int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
     }
     nrows = G_window_rows();
     ncols = G_window_cols();
-    buffer = G_allocate_cell_buf();
-    dbuffer = (double *)G_malloc(ncols * sizeof(double));
+    dbuffer = G_allocate_d_raster_buf();
     for (row = 0; row < nrows; row++) {
-	if (G_get_c_raster_row(map_fd, buffer, row) < 0) {
-	    G_free(buffer);
+	if (G_get_d_raster_row(map_fd, dbuffer, row) < 0) {
 	    G_free(dbuffer);
 	    G_close_cell(map_fd);
 	    sprintf(msg, "%s(): unable to read file [%s] in [%s], %d %d",
@@ -36,11 +33,7 @@ int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
 	    G_warning(msg);
 	    return -2;
 	}
-	for (col = ncols - 1; col >= 0; col--) {
-	    dbuffer[col] = (double)buffer[col];
-	}
-	if (segment_put_row(&(dseg->seg), (CELL *) dbuffer, row) < 0) {
-	    G_free(buffer);
+	if (segment_put_row(&(dseg->seg), (DCELL *) dbuffer, row) < 0) {
 	    G_free(dbuffer);
 	    G_close_cell(map_fd);
 	    sprintf(msg, "%s(): unable to segment put row for [%s] in [%s]",
@@ -51,7 +44,6 @@ int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
     }
 
     G_close_cell(map_fd);
-    G_free(buffer);
     G_free(dbuffer);
 
     dseg->name = G_store(map_name);

raster/r.watershed/seg/dseg_write.c

@@ -7,27 +7,21 @@ static char *me = "dseg_write_cell";
 int dseg_write_cellfile(DSEG * dseg, char *map_name)
 {
     int map_fd;
-    int row, col, nrows, ncols;
-    CELL *buffer;
+    int row, nrows, ncols;
     double *dbuffer;
 
-    map_fd = G_open_cell_new(map_name);
+    map_fd = G_open_raster_new(map_name, DCELL_TYPE);
     if (map_fd < 0) {
 	G_warning("%s(): unable to open new map layer [%s]", me, map_name);
 	return -1;
     }
     nrows = G_window_rows();
     ncols = G_window_cols();
-    buffer = G_allocate_cell_buf();
-    dbuffer = (double *)G_malloc(ncols * sizeof(double));
+    dbuffer = G_allocate_d_raster_buf();
     segment_flush(&(dseg->seg));
     for (row = 0; row < nrows; row++) {
-	segment_get_row(&(dseg->seg), (CELL *) dbuffer, row);
-	for (col = ncols - 1; col >= 0; col--) {
-	    buffer[col] = (CELL) (dbuffer[col] + 0.5);
-	}
-	if (G_put_raster_row(map_fd, buffer, CELL_TYPE) < 0) {
-	    G_free(buffer);
+	segment_get_row(&(dseg->seg), (DCELL *) dbuffer, row);
+	if (G_put_raster_row(map_fd, dbuffer, DCELL_TYPE) < 0) {
 	    G_free(dbuffer);
 	    G_unopen_cell(map_fd);
 	    G_warning("%s(): unable to write new map layer [%s], row %d",
@@ -35,7 +29,6 @@ int dseg_write_cellfile(DSEG * dseg, char *map_name)
 	    return -2;
 	}
     }
-    G_free(buffer);
     G_free(dbuffer);
     G_close_cell(map_fd);
     return 0;

raster/r.watershed/seg/init_vars.c

@@ -15,7 +15,7 @@ int init_vars(int argc, char *argv[])
 
     /* int page_block, num_cseg; */
     int max_bytes;
-    CELL *buf, alt_value, asp_value, worked_value;
+    CELL *buf, alt_value, asp_value, worked_value, block_value;
     DCELL wat_value;
     char MASK_flag;
 
@@ -131,8 +131,7 @@ int init_vars(int argc, char *argv[])
 
     /* segment parameters: one size fits all. Fine tune? */
     /* Segment rows and cols: 200 */
-    /* 1 segment open for all rasters: 2.86 MB */
-    /* num_open_segs = segs_mb / 2.86 */
+    /* 1 segment open for all rasters: 1.34 MB */
 
     seg_rows = SROW;
     seg_cols = SCOL;
@@ -142,7 +141,7 @@ int init_vars(int argc, char *argv[])
 	G_warning(_("Maximum memory to be used was smaller than 3 MB, set to default = 300 MB."));
     }
 
-    num_open_segs = segs_mb / 2.86;
+    num_open_segs = segs_mb / 1.34;
 
     G_debug(1, "segs MB: %.0f", segs_mb);
     G_debug(1, "region rows: %d", nrows);
@@ -166,24 +165,92 @@ int init_vars(int argc, char *argv[])
     G_debug(1, "  open segments after adjusting:\t%d", num_open_segs);
 
     cseg_open(&alt, seg_rows, seg_cols, num_open_segs);
-    cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
     cseg_read_cell(&alt, ele_name, "");
-    cseg_read_cell(&r_h, ele_name, "");
+    if (er_flag) {
+	cseg_open(&r_h, seg_rows, seg_cols, num_open_segs);
+	cseg_read_cell(&r_h, ele_name, "");
+    }
+    
+    /* NULL cells in input elevation map */
+    bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
+    bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
+    G_debug(1, "Checking for masked and NULL cells in input elevation <%s>.", ele_name);
+    MASK_flag = 0;
+    do_points = nrows * ncols;
+    fd = G_open_cell_old(ele_name, "");
+    if (fd < 0) {
+	G_fatal_error(_("unable to open elevation map layer"));
+    }
+    buf = G_allocate_cell_buf();
+    for (r = 0; r < nrows; r++) {
+	G_get_c_raster_row(fd, buf, r);
+	for (c = 0; c < ncols; c++) {
+	    /* check for masked and NULL cells */
+	    alt_value = buf[c];
+	    if (G_is_c_null_value(&alt_value)) {
+		bseg_put(&worked, &one, r, c);
+		bseg_put(&in_list, &one, r, c);
+		do_points--;
+	    }
+	}
+    }
+    G_close_cell(fd);
+    G_free(buf);
+    if (do_points < nrows * ncols)
+	MASK_flag = 1;
+    
+    /* initial flow accumulation */
     dseg_open(&wat, seg_rows, seg_cols, num_open_segs);
-
     if (run_flag) {
 	dseg_read_cell(&wat, run_name, "");
+	if (MASK_flag) {
+	    for (r = 0; r < nrows; r++) {
+		for (c = 0; c < ncols; c++) {
+		    bseg_get(&worked, &worked_value, r, c);
+		    if (worked_value)
+			dseg_put(&wat, &d_zero, r, c);
+		}
+	    }
+	}
     }
     else {
 	for (r = 0; r < nrows; r++) {
 	    for (c = 0; c < ncols; c++)
-		if (-1 == dseg_put(&wat, &d_one, r, c))
-		    exit(EXIT_FAILURE);
+		if (MASK_flag) {
+		    bseg_get(&worked, &worked_value, r, c);
+		    if (worked_value)
+			dseg_put(&wat, &d_zero, r, c);
+		    else
+			dseg_put(&wat, &d_one, r, c);
+		}
+		else {
+		    if (-1 == dseg_put(&wat, &d_one, r, c))
+			exit(EXIT_FAILURE);
+		}
 	}
     }
     cseg_open(&asp, seg_rows, seg_cols, num_open_segs);
+    /* depression: drainage direction will be set to zero later */
     if (pit_flag) {
-	cseg_read_cell(&asp, pit_name, "");
+	fd = G_open_cell_old(pit_name, "");
+	if (fd < 0) {
+	    G_fatal_error(_("unable to open depression map layer"));
+	}
+	buf = G_allocate_cell_buf();
+	for (r = 0; r < nrows; r++) {
+	    G_get_c_raster_row(fd, buf, r);
+	    for (c = 0; c < ncols; c++) {
+		asp_value = buf[c];
+		if (!G_is_c_null_value(&asp_value)) {
+		    cseg_put(&asp, &one, r, c);
+		}
+		else {
+		    cseg_put(&asp, &zero, r, c);
+		}
+	    }
+	}
+	G_close_cell(fd);
+	G_free(buf);
     }
     else {
 	for (r = 0; r < nrows; r++) {
@@ -194,7 +261,25 @@ int init_vars(int argc, char *argv[])
     }
     bseg_open(&swale, seg_rows, seg_cols, num_open_segs);
     if (ob_flag) {
-	bseg_read_cell(&swale, ob_name, "");
+	fd = G_open_cell_old(ob_name, "");
+	if (fd < 0) {
+	    G_fatal_error(_("unable to open blocking map layer"));
+	}
+	buf = G_allocate_cell_buf();
+	for (r = 0; r < nrows; r++) {
+	    G_get_c_raster_row(fd, buf, r);
+	    for (c = 0; c < ncols; c++) {
+		block_value = buf[c];
+		if (!G_is_c_null_value(&block_value)) {
+		    bseg_put(&swale, &one, r, c);
+		}
+		else {
+		    bseg_put(&swale, &zero, r, c);
+		}
+	    }
+	}
+	G_close_cell(fd);
+	G_free(buf);
     }
     else {
 	for (r = 0; r < nrows; r++) {
@@ -206,44 +291,26 @@ int init_vars(int argc, char *argv[])
 	dseg_open(&ril, 1, seg_rows * seg_cols, num_open_segs);
 	dseg_read_cell(&ril, ril_name, "");
     }
-    bseg_open(&in_list, seg_rows, seg_cols, num_open_segs);
-    bseg_open(&worked, seg_rows, seg_cols, num_open_segs);
-    MASK_flag = 0;
-    do_points = nrows * ncols;
-    if (NULL != G_find_file("cell", "MASK", G_mapset())) {
-	MASK_flag = 1;
-	if ((fd = G_open_cell_old("MASK", G_mapset())) < 0) {
-	    G_fatal_error(_("Unable to open MASK"));
-	}
-	else {
-	    buf = G_allocate_cell_buf();
-	    for (r = 0; r < nrows; r++) {
-		G_get_c_raster_row_nomask(fd, buf, r);
-		for (c = 0; c < ncols; c++) {
-		    if (!buf[c]) {
-			do_points--;
-			bseg_put(&worked, &one, r, c);
-			bseg_put(&in_list, &one, r, c);
-		    }
-		}
-	    }
-	    G_close_cell(fd);
-	    G_free(buf);
-	}
-    }
+    
     /* dseg_open(&slp, SROW, SCOL, num_open_segs); */
-    dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
+
+    /* RUSLE: LS and/or S factor */
+
+    if (er_flag) {
+	dseg_open(&s_l, seg_rows, seg_cols, num_open_segs);
+    }
     if (sg_flag)
 	dseg_open(&s_g, 1, seg_rows * seg_cols, num_open_segs);
     if (ls_flag)
 	dseg_open(&l_s, 1, seg_rows * seg_cols, num_open_segs);
-    seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols,
-	     num_open_segs, sizeof(POINT));
+
+    seg_open(&astar_pts, 1, do_points, 1, seg_rows * seg_cols * 2,
+	     num_open_segs / 2, sizeof(POINT));
 
     /* heap_index will track astar_pts in ternary min-heap */
     /* heap_index is one-based */
-    seg_open(&heap_index, 1, do_points + 1, 1, seg_rows * seg_cols,
-	     num_open_segs, sizeof(HEAP));
+    seg_open(&heap_index, 1, do_points + 1, 1, seg_cols * num_open_segs * seg_rows / 10,
+	     10, sizeof(HEAP));
 
     G_message(_("SECTION 1b (of %1d): Determining Offmap Flow."), tot_parts);
 
@@ -261,7 +328,8 @@ int init_vars(int argc, char *argv[])
 		    dseg_put(&wat, &d_zero, r, c);
 		}
 		else {
-		    dseg_put(&s_l, &half_res, r, c);
+		    if (er_flag)
+			dseg_put(&s_l, &half_res, r, c);
 		    cseg_get(&asp, &asp_value, r, c);
 		    if (r == 0 || c == 0 || r == nrows - 1 ||
 			c == ncols - 1 || asp_value != 0) {
@@ -270,7 +338,15 @@ int init_vars(int argc, char *argv[])
 			    wat_value = -wat_value;
 			    dseg_put(&wat, &wat_value, r, c);
 			}
-			if (r == 0)
+			/* set depression */
+			if (asp_value) {
+			    asp_value = 0;
+			    if (wat_value < 0) {
+				wat_value = -wat_value;
+				dseg_put(&wat, &wat_value, r, c);
+			    }
+			}
+			else if (r == 0)
 			    asp_value = -2;
 			else if (c == 0)
 			    asp_value = -4;
@@ -278,17 +354,15 @@ int init_vars(int argc, char *argv[])
 			    asp_value = -6;
 			else if (c == ncols - 1)
 			    asp_value = -8;
-			else
-			    asp_value = -1;
 			if (-1 == cseg_put(&asp, &asp_value, r, c))
 			    exit(EXIT_FAILURE);
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 		    }
 		    else if (!bseg_get(&worked, &worked_value, r - 1, c)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -2;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -300,7 +374,7 @@ int init_vars(int argc, char *argv[])
 		    else if (!bseg_get(&worked, &worked_value, r + 1, c)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -6;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -312,7 +386,7 @@ int init_vars(int argc, char *argv[])
 		    else if (!bseg_get(&worked, &worked_value, r, c - 1)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -4;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -324,7 +398,7 @@ int init_vars(int argc, char *argv[])
 		    else if (!bseg_get(&worked, &worked_value, r, c + 1)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -8;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -337,7 +411,7 @@ int init_vars(int argc, char *argv[])
 			     !bseg_get(&worked, &worked_value, r - 1, c - 1)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -3;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -350,7 +424,7 @@ int init_vars(int argc, char *argv[])
 			     !bseg_get(&worked, &worked_value, r - 1, c + 1)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -1;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -363,7 +437,7 @@ int init_vars(int argc, char *argv[])
 			     !bseg_get(&worked, &worked_value, r + 1, c - 1)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -5;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -376,7 +450,7 @@ int init_vars(int argc, char *argv[])
 			     !bseg_get(&worked, &worked_value, r + 1, c + 1)
 			     && worked_value != 0) {
 			cseg_get(&alt, &alt_value, r, c);
-			add_pt(r, c, -1, -1, alt_value, alt_value);
+			add_pt(r, c, alt_value, alt_value);
 			asp_value = -7;
 			cseg_put(&asp, &asp_value, r, c);
 			dseg_get(&wat, &wat_value, r, c);
@@ -398,7 +472,8 @@ int init_vars(int argc, char *argv[])
 	    G_percent(r, nrows, 2);
 	    for (c = 0; c < ncols; c++) {
 		bseg_put(&worked, &zero, r, c);
-		dseg_put(&s_l, &half_res, r, c);
+		if (er_flag)
+		    dseg_put(&s_l, &half_res, r, c);
 		cseg_get(&asp, &asp_value, r, c);
 		if (r == 0 || c == 0 || r == nrows - 1 ||
 		    c == ncols - 1 || asp_value != 0) {
@@ -408,7 +483,15 @@ int init_vars(int argc, char *argv[])
 			if (-1 == dseg_put(&wat, &wat_value, r, c))
 			    exit(EXIT_FAILURE);
 		    }
-		    if (r == 0)
+		    /* set depression */
+		    if (asp_value) {
+			asp_value = 0;
+			if (wat_value < 0) {
+			    wat_value = -wat_value;
+			    dseg_put(&wat, &wat_value, r, c);
+			}
+		    }
+		    else if (r == 0)
 			asp_value = -2;
 		    else if (c == 0)
 			asp_value = -4;
@@ -416,12 +499,10 @@ int init_vars(int argc, char *argv[])
 			asp_value = -6;
 		    else if (c == ncols - 1)
 			asp_value = -8;
-		    else
-			asp_value = -1;
 		    if (-1 == cseg_put(&asp, &asp_value, r, c))
 			exit(EXIT_FAILURE);
 		    cseg_get(&alt, &alt_value, r, c);
-		    add_pt(r, c, -1, -1, alt_value, alt_value);
+		    add_pt(r, c, alt_value, alt_value);
 		}
 		else {
 		    bseg_put(&in_list, &zero, r, c);

raster/r.watershed/seg/main.c

@@ -47,7 +47,6 @@ SHORT nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
 char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 char run_name[GNAME_MAX], ob_name[GNAME_MAX];
 char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
-
 const char *this_mapset;
 char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX],
     thr_name[8];

raster/r.watershed/seg/sg_factor.c

@@ -9,10 +9,9 @@ int sg_factor(void)
     CELL downer, low_elev, hih_elev;
     double height, length, S, sin_theta;
 
-    G_message(_("SECTION 4: Length Slope determination."));
+    G_message(_("SECTION 4: RUSLE LS and/or S factor determination."));
     for (r = nrows - 1; r >= 0; r--) {
-/* FIXME: G_percent params count backwards */
-	G_percent(r, nrows, 3);
+	G_percent(nrows - r, nrows, 3);
 	for (c = ncols - 1; c >= 0; c--) {
 	    cseg_get(&alt, &low_elev, r, c);
 	    cseg_get(&r_h, &hih_elev, r, c);
@@ -38,7 +37,7 @@ int sg_factor(void)
 	    }
 	}
     }
-    G_percent(r, nrows, 3);	/* finish it */
+    G_percent(nrows, nrows, 1);	/* finish it */
 
     return 0;
 }

raster/r.watershed/seg/slope_len.c

@@ -20,7 +20,7 @@ int slope_length(SHORT r, SHORT c, SHORT dr, SHORT dc)
 	    if (asp_value == 2 || asp_value == 6)
 		res = window.ns_res;
 	    else		/* asp_value == 4, 8, -2, -4, -6, or -8 */
-		res = diag;
+		res = diag;     /* how can res be diag with sides == 4??? */
 	}
 	else {			/* c == dc */