/***********************************
/* Program: ls.aml
/*
/* Function: Calculate a LS value for each field in the subwatershed
/* to use as input to ARDBSN file.
/*
/* History: Created 1999 by Jill Heller, the Advanced Resource Technology
/* Group, the University of Arizona, Tucson, AZ 85721.
/* ph: 520-621-3045. jheller@nexus.srnr.arizona.edu
/***********************************
/* Routines called: None.
/*
/* Requirements: Elevation GRID (DEM) and subwatershed coverage.
/*
/* Results: USLE length-slope (LS) factor map containing
/* values for each field of the subwatershed.
/*
/* Usage: &run ls
/*
/* Arguments: None
/**********************************
&IF %:PROGRAM% ne GRID &THEN &do
&IF %:PROGRAM% ne ARC &THEN quit
GRID
&end
&do sws &list <subws>
%sws%_grd = polygrid (%sws%, #, #, #, 30)
/* the above converted the sub-basin coverage to a grid with a cell
/* size = 30m.
/* To calculate area, the COUNT can be multiplied by 30*30 = 900m^2
/* for each cell in the COUNT.
/* The grid comes up with 60 sub-basins. This seems right for now.
gridclip dem6_usgs %sws%_elev cover %sws%
%sws%_slp = slope (%sws%_elev, percentrise)
/* now have a slope map for the entire area encompassing sub-basin CEDAR
/* next, will clip the slope map to the subwatershed coverage
/* Derive the average slope for each subwatershed
%sws%_av = zonalmean (%sws%_grd, %sws%_slp)
/* sws_av grid will serve as the "s" layer to the USLE
/*
gridclip dem6_dir fldir_%sws% cover %sws%
/* the fldir_usgs grid is the dem clipped to the watershed boundary
/* w/ usgs gage as outlet.
/*
fl_%sws% = flowlength (fldir_%sws%)
/*calcs flowlength from flowdir
/*
max_%sws% = zonalmax (%sws%_grd, fl_%sws%)
min_%sws% = zonalmin (%sws%_grd, fl_%sws%)
%sws%_lngth = max_%sws% - min_%sws%
/*cedr_length will serve as "L" factor
L_eng%sws% = 3.208 * %sws%_lngth
/* converts the sws_length grid to English units from metric (meters
/* to feet)
L_%sws% = sqrt (L_eng%sws% div 72.6)
/* creates the "L" factor grid to be used in the USLE equation.
ssqr_%sws% = sqr (%sws%_av)
s_%sws% = (.43 + (.3 * %sws%_av) + (.043 * ssqr_%sws%)) div 6.613
LS_%sws% = L_%sws% * (s_%sws% * (10000 div (10000 + sqr(s_%sws%))))
kill (!l_%sws%, ssqr_%sws%, s_%sws%, l_eng%sws%, %sws%_lngth, max_%sws%, min_%sws%!) all
/* mapex %sws%
/* gridpaint %sws%
/* CELLVALUE LS_%sws% *
/* Extract the length-slope (LS) factor values for each field by clicking in
/* each field and recording the given value.
&end
&return