r.watershed: +retention for flow distribution, credits to Andreas Gericke (IGB Berlin)

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

raster/r.watershed/front/main.c

@@ -60,6 +60,7 @@ int main(int argc, char *argv[])
     struct Option *opt16;
     struct Option *opt17;
     struct Option *opt18;
+    struct Option *opt19;
     struct Flag *flag_sfd;
     struct Flag *flag_flow;
     struct Flag *flag_seg;
@@ -116,6 +117,13 @@ int main(int argc, char *argv[])
     opt5->required = NO;
     opt5->guisection = _("Inputs");
 
+    opt19 = G_define_standard_option(G_OPT_R_INPUT);
+    opt19->key = "retention";
+    opt19->description =
+        _("Name of input raster map with percentages for flow accumulation.");
+    opt19->required = NO;
+    opt19->guisection = _("Inputs");
+
     opt6 = G_define_option();
     opt6->key = "threshold";
     opt6->description = _("Minimum size of exterior watershed basin");
@@ -306,6 +314,7 @@ int main(int argc, char *argv[])
     do_opt(opt3);
     do_opt(opt4);
     do_opt(opt5);
+    do_opt(opt19);
     do_opt(opt6);
     do_opt(opt7);
     do_opt(opt8);

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

@@ -73,6 +73,12 @@ represent the amount of overland flow each cell contributes to surface
 flow. If omitted, a value of one (1) is assumed.
 
 <p>
+Raster <b>retention</b> map specifies correction factors per cell for 
+flow distribution. This map indicates the percentage of overland flow 
+units leaving each cell. Values should be between zero and 100. If 
+omitted, a value of 100 is assumed.
+
+<p>
 Input Raster map or value containing the percent of disturbed land
 (i.e., croplands, and construction sites) where the raster or input
 value of 17 equals 17%.  If no map or value is
@@ -567,6 +573,8 @@ Charles Ehlschlaeger, U.S. Army Construction Engineering Research Laboratory
 <br>
 Faster sorting algorithm and MFD support:
 Markus Metz &lt;markus.metz.giswork at gmail.com&gt;
+<br>
+Retention for flow distribution by Andreas Gericke (IGB Berlin)
 
 <p>
 <i>Last changed: $Date$</i>

raster/r.watershed/ram/Gwater.h

@@ -56,10 +56,11 @@ extern OC_STACK *ocs;
 extern int ocs_alloced;
 extern FLAG *worked, *in_list, *s_b, *swale, *flat_done;
 extern RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
-extern RAMSEG r_h_seg, dep_seg;
+extern RAMSEG r_h_seg, dep_seg, rtn_seg;
 extern RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
 extern int *astar_pts;
 extern CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+extern char *rtn;
 extern DCELL *wat, *sca, *tanb;
 extern int ril_fd;
 extern double *s_l, *s_g, *l_s;
@@ -72,13 +73,13 @@ extern int nextdr[8];
 extern int nextdc[8];
 extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
-extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
+extern char ril_name[GNAME_MAX], rtn_name[GNAME_MAX], dep_name[GNAME_MAX];
 extern const char *this_mapset;
 extern char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
 extern char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[GNAME_MAX];
 extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX], spi_name[GNAME_MAX];
 extern char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
-extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag;
+extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag, rtn_flag;
 extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, tci_flag, spi_flag, atanb_flag;
 extern char bas_flag, seg_flag, haf_flag, er_flag;
 extern char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;

raster/r.watershed/ram/close_maps.c

@@ -21,6 +21,9 @@ int close_maps(void)
     if (wat_flag || ls_flag || sl_flag || sg_flag || atanb_flag)
 	dbuf = Rast_allocate_d_buf();
     G_free(alt);
+    if (rtn_flag)
+	G_free(rtn);
+
     if (ls_flag || sg_flag)
 	G_free(r_h);
 

raster/r.watershed/ram/do_cum.c

@@ -140,6 +140,10 @@ int do_cum(void)
 	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {	/* if ((dr = astar_pts[killer].downr) > -1) { */
 	    down_index = SEG_INDEX(wat_seg, dr, dc);
 	    value = wat[this_index];
+            /* apply retention to adjust flow accumulation */
+            if (rtn_flag)
+                value *= rtn[this_index] / 100.0;
+
 	    if (fabs(value) >= threshold)
 		FLAG_SET(swale, r, c);
 	    valued = wat[down_index];
@@ -197,7 +201,7 @@ int do_cum(void)
 	    /* topographic wetness index ln(a / tan(beta)) and
 	     * stream power index a * tan(beta) */
 	    if (atanb_flag) {
-		sca[this_index] = fabs(wat[this_index]) *
+		sca[this_index] = fabs(value) *
 		    (cell_size / contour[np_side]);
 		tanb[this_index] = get_slope_tci(alt[this_index],
 						 alt[down_index],
@@ -321,6 +325,9 @@ int do_cum_mfd(void)
 	    dr = dc = -1;
 	if (dr >= 0 && dr < nrows && dc >= 0 && dc < ncols) {	/* if ((dr = astar_pts[killer].downr) > -1) { */
 	    value = wat[this_index];
+            /* apply retention to adjust flow accumulation */
+            if (rtn_flag)
+                value *= rtn[this_index] / 100.0;
 	    down_index = SEG_INDEX(wat_seg, dr, dc);
 
 	    /* get weights */
@@ -489,7 +496,7 @@ int do_cum_mfd(void)
 	    /* topographic wetness index ln(a / tan(beta)) and
 	     * stream power index a * tan(beta) */
 	    if (atanb_flag) {
-		sca[this_index] = fabs(wat[this_index]) *
+		sca[this_index] = fabs(value) *
 		    (cell_size / sum_contour);
 		tanb[this_index] = tci_div;
 	    }

raster/r.watershed/ram/init_vars.c

@@ -20,7 +20,7 @@ int init_vars(int argc, char *argv[])
 
     G_gisinit(argv[0]);
     /* input */
-    ele_flag = pit_flag = run_flag = ril_flag = 0;
+    ele_flag = pit_flag = run_flag = ril_flag = rtn_flag = 0;
     /* output */
     wat_flag = asp_flag = tci_flag = spi_flag = atanb_flag = 0;
     bas_flag = seg_flag = haf_flag = 0;
@@ -67,6 +67,8 @@ int init_vars(int argc, char *argv[])
 	    haf_flag++;
 	else if (sscanf(argv[r], "flow=%s", run_name) == 1)
 	    run_flag++;
+	else if (sscanf(argv[r], "retention=%s", rtn_name) == 1)
+	    rtn_flag++;
 	else if (sscanf(argv[r], "ar=%s", arm_name) == 1)
 	    arm_flag++;
 	/* slope length
@@ -266,6 +268,42 @@ int init_vars(int argc, char *argv[])
 	G_free(dbuf);
     }
 
+    /* read retention map to adjust flow distribution (AG) */
+    rtn = NULL;
+    if (rtn_flag) {
+	rtn = (char *) G_malloc(sizeof(char) *
+                           size_array(&rtn_seg, nrows, ncols));
+        buf = Rast_allocate_c_buf();
+        fd = Rast_open_old(rtn_name, "");
+        for (r = 0; r < nrows; r++) {
+            Rast_get_c_row(fd, buf, r);
+            for (c = 0; c < ncols; c++) {
+                if (MASK_flag) {
+                    block_value = FLAG_GET(worked, r, c);
+                    if (!block_value) {
+                        block_value = buf[c];
+                    }
+		    else
+			block_value = 100;
+                }
+                else
+                    block_value = buf[c];
+
+                if (!Rast_is_c_null_value(&block_value)) {
+		    if (block_value < 0)
+			block_value = 0;
+		    if (block_value > 100)
+			block_value = 100;
+                    rtn[SEG_INDEX(rtn_seg, r, c)] = block_value;
+		}
+                else
+                    rtn[SEG_INDEX(rtn_seg, r, c)] = 100;
+            }
+        }
+        Rast_close(fd);
+        G_free(buf);
+    }
+
     /* overland blocking map; this is also creating streams... */
     if (ob_flag) {
 	buf = Rast_allocate_c_buf();

raster/r.watershed/ram/main.c

@@ -35,10 +35,11 @@ OC_STACK *ocs;
 int ocs_alloced;
 FLAG *worked, *in_list, *s_b, *swale, *flat_done;
 RAMSEG dis_seg, alt_seg, wat_seg, asp_seg, bas_seg, haf_seg;
-RAMSEG r_h_seg, dep_seg;
+RAMSEG r_h_seg, dep_seg, rtn_seg;
 RAMSEG slp_seg, s_l_seg, s_g_seg, l_s_seg;
 int *astar_pts;
 CELL *dis, *alt, *asp, *bas, *haf, *r_h, *dep;
+char *rtn;
 DCELL *wat, *sca, *tanb;
 int ril_fd;
 double *s_l, *s_g, *l_s;
@@ -52,7 +53,7 @@ int 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];
+char ril_name[GNAME_MAX], rtn_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];
@@ -62,7 +63,7 @@ char wat_name[GNAME_MAX], asp_name[GNAME_MAX], tci_name[GNAME_MAX],
     spi_name[GNAME_MAX];
 char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
 char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, flat_flag;
-char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, rtn_flag;
 char bas_flag, seg_flag, haf_flag, er_flag, tci_flag, spi_flag, atanb_flag;
 char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
 FILE *fp;

raster/r.watershed/ram/usage.c

@@ -10,11 +10,12 @@ void usage(char *me)
 		    "threshold=swale_threshold [flow=overland_flow_map] "
 		    "[drainage=drain_direction_map] [depression=depression_map] "
 		    "[accumulation=accumulation_map] [basin=watershed_basin_map] "
-		    "[stream=stream_segment_map]\n\n"
+		    "[retention=retention_map] [stream=stream_segment_map]\n\n"
 		    "USAGE for slope length determination:\n%s [-4] "
 		    "elevation=elevation_map threshold=swale_threshold "
 		    "[drainage=drain_direction_map] [depression=depression_map] "
-		    "[accumulation=accumulation_map] [max_slope_length=max_slope_length] "
+		    "[accumulation=accumulation_map] [retention=retention_map] "
+		    "[max_slope_length=max_slope_length] "
 		    "[blocking=overland_blocking_map] [slope_steepness=slope_steepness_map] "
 		    "length_slope=length_slope_map [disturbed_land=rill_erosion_map] "
 		    "[slope_deposition=slope_deposition value or map]"
@@ -23,5 +24,5 @@ void usage(char *me)
 		    "[drainage=drain_direction_map] [depression=depression_map] "
 		    "[accumulation=accumulation_map] [basin=watershed_basin_map] "
 		    "[stream=stream_segment_map] [half_basin=half_basin_map] "
-		    "ar=ARMSED_file_name\n\n"), me, me, me);
+		    "[retention=retention_map] ar=ARMSED_file_name\n\n"), me, me, me);
 }

raster/r.watershed/seg/Gwater.h

@@ -93,7 +93,7 @@ extern int ocs_alloced;
 extern double half_res, diag, max_length, dep_slope;
 extern int bas_thres, tot_parts;
 extern SSEG astar_pts;
-extern BSEG s_b;
+extern BSEG s_b, rtn;
 extern CSEG dis, bas, haf, r_h, dep;
 extern SSEG watalt, aspflag;
 extern DSEG slp, s_l, s_g, l_s, ril;
@@ -108,7 +108,7 @@ extern int nextdr[8];
 extern int nextdc[8];
 extern char ele_name[GNAME_MAX], pit_name[GNAME_MAX];
 extern char run_name[GNAME_MAX], ob_name[GNAME_MAX];
-extern char ril_name[GNAME_MAX], dep_name[GNAME_MAX];
+extern char ril_name[GNAME_MAX], rtn_name[GNAME_MAX], dep_name[GNAME_MAX];
 
 extern const char *this_mapset;
 extern char seg_name[GNAME_MAX], bas_name[GNAME_MAX], haf_name[GNAME_MAX], thr_name[8];
@@ -116,7 +116,7 @@ extern char ls_name[GNAME_MAX], st_name[GNAME_MAX], sl_name[GNAME_MAX], sg_name[
 extern char wat_name[GNAME_MAX], asp_name[GNAME_MAX];
 extern char tci_name[GNAME_MAX], spi_name[GNAME_MAX];
 extern char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
-extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
+extern char ele_flag, pit_flag, run_flag, dis_flag, ob_flag, rtn_flag;
 extern char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, tci_flag, spi_flag, atanb_flag;
 extern char bas_flag, seg_flag, haf_flag, er_flag;
 extern char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;

raster/r.watershed/seg/bseg_read.c

@@ -5,8 +5,8 @@
 
 int bseg_read_cell(BSEG * bseg, char *map_name, char *mapset)
 {
-    int row, nrows;
-    int col, ncols;
+    int row, rows;
+    int col, cols;
     int map_fd;
     CELL *buffer;
     char cbuf;
@@ -15,12 +15,12 @@ int bseg_read_cell(BSEG * bseg, char *map_name, char *mapset)
     bseg->mapset = NULL;
 
     map_fd = Rast_open_old(map_name, mapset);
-    nrows = Rast_window_rows();
-    ncols = Rast_window_cols();
+    rows = Rast_window_rows();
+    cols = Rast_window_cols();
     buffer = Rast_allocate_c_buf();
-    for (row = 0; row < nrows; row++) {
+    for (row = 0; row < rows; row++) {
 	Rast_get_c_row(map_fd, buffer, row);
-	for (col = ncols; col >= 0; col--) {
+	for (col = cols; col >= 0; col--) {
 	    cbuf = (char)buffer[col];
 	    bseg_put(bseg, &cbuf, row, col);
 	}

raster/r.watershed/seg/bseg_write.c

@@ -5,24 +5,24 @@
 int bseg_write_cellfile(BSEG * bseg, char *map_name)
 {
     int map_fd;
-    int row, nrows;
-    int col, ncols;
+    int row, rows;
+    int col, cols;
     CELL *buffer;
     char value;
 
     map_fd = Rast_open_c_new(map_name);
-    nrows = Rast_window_rows();
-    ncols = Rast_window_cols();
+    rows = Rast_window_rows();
+    cols = Rast_window_cols();
     buffer = Rast_allocate_c_buf();
-    for (row = 0; row < nrows; row++) {
-	G_percent(row, nrows, 1);
-	for (col = 0; col < ncols; col++) {
+    for (row = 0; row < rows; row++) {
+	G_percent(row, rows, 1);
+	for (col = 0; col < cols; col++) {
 	    bseg_get(bseg, &value, row, col);
 	    buffer[col] = value;
 	}
 	Rast_put_row(map_fd, buffer, CELL_TYPE);
     }
-    G_percent(row, nrows, 1);	/* finish it */
+    G_percent(row, rows, 1);	/* finish it */
     G_free(buffer);
     Rast_close(map_fd);
     return 0;

raster/r.watershed/seg/close_maps.c

@@ -16,6 +16,9 @@ int close_maps(void)
     WAT_ALT *wabuf;
     ASP_FLAG *afbuf;
 
+    if (rtn_flag)
+        bseg_close(&rtn);
+
     if (wat_flag) {
 	G_message(_("Closing accumulation map"));
 	sum = sum_sqr = stddev = 0.0;

raster/r.watershed/seg/cseg_read.c

@@ -7,7 +7,7 @@ static char *me = "cseg_read_cell";
 
 int cseg_read_cell(CSEG * cseg, char *map_name, char *mapset)
 {
-    GW_LARGE_INT row, nrows;
+    GW_LARGE_INT row, rows;
     int map_fd;
     CELL *buffer;
 
@@ -15,9 +15,9 @@ int cseg_read_cell(CSEG * cseg, char *map_name, char *mapset)
     cseg->mapset = NULL;
 
     map_fd = Rast_open_old(map_name, mapset);
-    nrows = Rast_window_rows();
+    rows = Rast_window_rows();
     buffer = Rast_allocate_c_buf();
-    for (row = 0; row < nrows; row++) {
+    for (row = 0; row < rows; row++) {
 	Rast_get_c_row(map_fd, buffer, row);
 	if (Segment_put_row(&(cseg->seg), buffer, row) < 0) {
 	    G_free(buffer);

raster/r.watershed/seg/cseg_write.c

@@ -6,19 +6,19 @@
 int cseg_write_cellfile(CSEG * cseg, char *map_name)
 {
     int map_fd;
-    GW_LARGE_INT row, nrows;
+    GW_LARGE_INT row, rows;
     CELL *buffer;
 
     map_fd = Rast_open_c_new(map_name);
-    nrows = Rast_window_rows();
+    rows = Rast_window_rows();
     buffer = Rast_allocate_c_buf();
     Segment_flush(&(cseg->seg));
-    for (row = 0; row < nrows; row++) {
-	G_percent(row, nrows, 1);
+    for (row = 0; row < rows; row++) {
+	G_percent(row, rows, 1);
 	Segment_get_row(&(cseg->seg), buffer, row);
 	Rast_put_row(map_fd, buffer, CELL_TYPE);
     }
-    G_percent(row, nrows, 1);	/* finish it */
+    G_percent(row, rows, 1);	/* finish it */
     G_free(buffer);
     Rast_close(map_fd);
     return 0;

raster/r.watershed/seg/do_cum.c

@@ -113,6 +113,7 @@ int do_cum(void)
     int asp_r[9] = { 0, -1, -1, -1, 0, 1, 1, 1, 0 };
     int asp_c[9] = { 0, 1, 0, -1, -1, -1, 0, 1, 1 };
     WAT_ALT wa, wadown;
+    char rtn_value;
     ASP_FLAG af, afdown;
     A_TANB sca_tanb;
     double *dist_to_nbr, *contour;
@@ -176,6 +177,10 @@ int do_cum(void)
 
 	    seg_get(&watalt, (char *)&wa, r, c);
 	    value = wa.wat;
+	    if (rtn_flag) {
+		bseg_get(&rtn, (char *)&rtn_value, dr, dc);
+		value *= rtn_value / 100.0;
+	    }
 	    is_swale = FLAG_GET(af.flag, SWALEFLAG);
 	    if (fabs(value) >= threshold && !is_swale) {
 		is_swale = 1;
@@ -201,7 +206,7 @@ int do_cum(void)
 	    /* topographic wetness index ln(a / tan(beta)) and
 	     * stream power index a * tan(beta) */
 	    if (atanb_flag) {
-		sca_tanb.sca = fabs(wa.wat) * (cell_size / contour[np_side]);
+		sca_tanb.sca = fabs(value) * (cell_size / contour[np_side]);
 
 		sca_tanb.tanb = get_slope_tci(wa.ele, wadown.ele,
 					      dist_to_nbr[np_side]);
@@ -276,6 +281,7 @@ int do_cum_mfd(void)
     double sum_contour, cell_size;
     POINT point;
     WAT_ALT wa;
+    char rtn_value;
     ASP_FLAG af, afdown;
     A_TANB sca_tanb;
     GW_LARGE_INT killer;
@@ -351,6 +357,10 @@ int do_cum_mfd(void)
 
 	    seg_get(&watalt, (char *)&wa, r, c);
 	    value = wa.wat;
+	    if (rtn_flag) {
+		bseg_get(&rtn, (char *)&rtn_value, dr, dc);
+		value *= rtn_value / 100.0;
+	    }
 
 	    /* get weights */
 	    max_weight = 0;

raster/r.watershed/seg/dseg_read.c

@@ -7,7 +7,7 @@ static char *me = "dseg_read_cell";
 
 int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
 {
-    GW_LARGE_INT row, nrows, ncols;
+    GW_LARGE_INT row, rows;
     int map_fd;
     double *dbuffer;
 
@@ -15,10 +15,9 @@ int dseg_read_cell(DSEG * dseg, char *map_name, char *mapset)
     dseg->mapset = NULL;
 
     map_fd = Rast_open_old(map_name, mapset);
-    nrows = Rast_window_rows();
-    ncols = Rast_window_cols();
+    rows = Rast_window_rows();
     dbuffer = Rast_allocate_d_buf();
-    for (row = 0; row < nrows; row++) {
+    for (row = 0; row < rows; row++) {
 	Rast_get_d_row(map_fd, dbuffer, row);
 	if (Segment_put_row(&(dseg->seg), (DCELL *) dbuffer, row) < 0) {
 	    G_free(dbuffer);

raster/r.watershed/seg/dseg_write.c

@@ -6,20 +6,19 @@
 int dseg_write_cellfile(DSEG * dseg, char *map_name)
 {
     int map_fd;
-    GW_LARGE_INT row, nrows, ncols;
+    GW_LARGE_INT row, rows;
     double *dbuffer;
 
     map_fd = Rast_open_new(map_name, DCELL_TYPE);
-    nrows = Rast_window_rows();
-    ncols = Rast_window_cols();
+    rows = Rast_window_rows();
     dbuffer = Rast_allocate_d_buf();
     Segment_flush(&(dseg->seg));
-    for (row = 0; row < nrows; row++) {
-	G_percent(row, nrows, 1);
+    for (row = 0; row < rows; row++) {
+	G_percent(row, rows, 1);
 	Segment_get_row(&(dseg->seg), (DCELL *) dbuffer, row);
 	Rast_put_row(map_fd, dbuffer, DCELL_TYPE);
     }
-    G_percent(row, nrows, 1);	/* finish it */
+    G_percent(row, rows, 1);	/* finish it */
     G_free(dbuffer);
     Rast_close(map_fd);
     return 0;

raster/r.watershed/seg/init_vars.c

@@ -32,7 +32,7 @@ int init_vars(int argc, char *argv[])
 
     G_gisinit(argv[0]);
     /* input */
-    ele_flag = pit_flag = run_flag = ril_flag = 0;
+    ele_flag = pit_flag = run_flag = ril_flag  = rtn_flag = 0;
     /* output */
     wat_flag = asp_flag = tci_flag = spi_flag = atanb_flag = 0;
     bas_flag = seg_flag = haf_flag = 0;
@@ -78,6 +78,8 @@ int init_vars(int argc, char *argv[])
 	    haf_flag++;
 	else if (sscanf(argv[r], "flow=%s", run_name) == 1)
 	    run_flag++;
+	else if (sscanf(argv[r], "retention=%s", rtn_name) == 1)
+	    rtn_flag++;
 	else if (sscanf(argv[r], "ar=%s", arm_name) == 1)
 	    arm_flag++;
 	/* slope length
@@ -261,6 +263,10 @@ int init_vars(int argc, char *argv[])
 	cseg_read_cell(&r_h, ele_name, "");
     }
 
+    if (rtn_flag) {
+	bseg_open(&rtn, seg_rows, seg_cols, num_open_segs);
+    }
+
     /* read elevation input and mark NULL/masked cells */
 
     /* scattered access: alt, watalt, bitflags, asp */
@@ -398,6 +404,34 @@ int init_vars(int argc, char *argv[])
 
     MASK_flag = (do_points < nrows * ncols);
 
+    /* read retention map to adjust flow distribution (AG) */
+    if (rtn_flag) {
+	char rtn_value;
+
+	fd = Rast_open_old(rtn_name, "");
+	buf = Rast_allocate_c_buf();
+	for (r = 0; r < nrows; r++) {
+	    G_percent(r, nrows, 1);
+	    Rast_get_c_row(fd, buf, r);
+	    for (c = 0; c < ncols; c++) {
+		block_value = buf[c];
+		if (Rast_is_c_null_value(&block_value))
+		    block_value = 100;
+		else {
+		    if (block_value < 0)
+			block_value = 0;
+		    if (block_value > 100)
+			block_value = 100;
+		}
+		rtn_value = block_value;
+		bseg_put(&rtn, &rtn_value, r, c);
+	    }
+	}
+	G_percent(nrows, nrows, 1);	/* finish it */
+	Rast_close(fd);
+	G_free(buf);
+    }
+
     /* do RUSLE */
     if (er_flag) {
 	if (ob_flag) {

raster/r.watershed/seg/main.c

@@ -36,7 +36,7 @@ int ocs_alloced;
 double half_res, diag, max_length, dep_slope;
 int bas_thres, tot_parts;
 SSEG astar_pts;
-BSEG s_b;
+BSEG s_b, rtn;
 CSEG dis, alt, bas, haf, r_h, dep;
 SSEG watalt, aspflag;
 DSEG slp, s_l, s_g, l_s, ril;
@@ -52,7 +52,7 @@ int 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];
+char ril_name[GNAME_MAX], rtn_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];
@@ -62,7 +62,7 @@ char wat_name[GNAME_MAX], asp_name[GNAME_MAX];
 char tci_name[GNAME_MAX], spi_name[GNAME_MAX];
 char arm_name[GNAME_MAX], dis_name[GNAME_MAX];
 char ele_flag, pit_flag, run_flag, dis_flag, ob_flag;
-char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag;
+char wat_flag, asp_flag, arm_flag, ril_flag, dep_flag, rtn_flag;
 char bas_flag, seg_flag, haf_flag, er_flag, tci_flag, spi_flag, atanb_flag;
 char st_flag, sb_flag, sg_flag, sl_flag, ls_flag;
 FILE *fp;

raster/r.watershed/seg/sseg_open.c

@@ -5,7 +5,7 @@
 #include "Gwater.h"
 
 int
-seg_open(SSEG * sseg, GW_LARGE_INT nrows, GW_LARGE_INT ncols, int row_in_seg,
+seg_open(SSEG * sseg, GW_LARGE_INT rows, GW_LARGE_INT cols, int row_in_seg,
 	 int col_in_seg, int nsegs_in_memory, int size_struct)
 {
     char *filename;
@@ -20,7 +20,7 @@ seg_open(SSEG * sseg, GW_LARGE_INT nrows, GW_LARGE_INT ncols, int row_in_seg,
 	G_warning("seg_open(): unable to create segment file");
 	return -2;
     }
-    if (0 > (errflag = Segment_format(fd, nrows, ncols,
+    if (0 > (errflag = Segment_format(fd, rows, cols,
 				      row_in_seg, col_in_seg, size_struct))) {
 	close(fd);
 	unlink(filename);