&args startgrid outgrid dirgrid elvgrd
/* Based on /puppet/esri/amls version, but looks at grid units
/* Harvey Greenberg hgreen@u.washington.edu
/* This runs very slowly
/* This takes a file grid of calculated lanslides, a direction grid, and
/* the flowdirection grid, and follows downstream runout paths.
/* Argument 1: grid of some unstable cells (mostly NODATA) All values must be 100.
/* Argument 2: Output grid: it will be 100/NODATA.  Existing grid will be deleted.
/* Argument 3: ESRI flowdirection grid
/* Argument 4 elevation grid (default scale factor = 1000)
/*
&s thresh_deg = 4 /* debris flows until the slope is 4 degrees
&s zscale 1000   /* If vertical units equal horizontal units, set this to 1

&describe %elvgrd%
&ty Using an angle of %thresh_deg% degrees and a Zscale of %zscale% of a %grd$dx%-unit grid
&setvar thresh  [calc [sin [angrad %thresh_deg%]] * %zscale% * %grd$dx%]
&setvar threshd [calc %thresh% * 1.414]
&ty Using internal threshholds of %thresh% and %threshd%
&setvar tempgrid t[after [date -ftag] .]
&type Our scratch file is %tempgrid%
&if [exists %tempgrid% -grid] &then; kill %tempgrid%
&if [exists %outgrid% -grid] &then; kill %outgrid%
copy %startgrid% %tempgrid%
&s oldmean = 0.9
&s its 1
/* big loop: until all holes are filled
&label bigloop
/* "else if"s do not work below.  Con statements fail in the same way.
  docell
    if(isnull(%tempgrid%)) {
      if((%tempgrid%(-1,0) > 0)  &  (%dirgrid%(-1,0) == 1) &~
      ((%elvgrd%(-1,0) - %elvgrd%) gt %thresh%)) %outgrid% = 1
      /* "else if" does not behave well
      if((%tempgrid%(-1,-1) > 0)  &  (%dirgrid%(-1,-1) ==   2) &~
         ((%elvgrd%(-1,-1) - %elvgrd%) gt %thresh%)) %outgrid% = 2
      if((%tempgrid%(0,-1) > 0)  &  (%dirgrid%( 0,-1) ==   4) &~
         ((%elvgrd%(0,-1) - %elvgrd%) gt %thresh%)) %outgrid% = 3
      if((%tempgrid%(1,-1) > 0)  &  (%dirgrid%( 1,-1) ==   8) &~
         ((%elvgrd%(1,-1) - %elvgrd%) gt %thresh%)) %outgrid% = 4
      if((%tempgrid%( 1, 0) > 0)  &  (%dirgrid%( 1, 0) ==  16) &~
         ((%elvgrd%(1,0) - %elvgrd%) gt %thresh%)) %outgrid% = 5
      if((%tempgrid%( 1, 1) > 0)  &  (%dirgrid%( 1, 1) ==  32) &~
         ((%elvgrd%( 1,1) - %elvgrd%) gt %thresh%)) %outgrid% = 6
      if((%tempgrid%( 0, 1) > 0)  &  (%dirgrid%( 0, 1) ==  64) &~
         ((%elvgrd%( 0,1) - %elvgrd%) gt %thresh%)) %outgrid% = 7
      if((%tempgrid%(-1, 1) > 0)  &  (%dirgrid%(-1, 1) == 128) &~
         ((%elvgrd%(-1,1) - %elvgrd%) gt %thresh%)) %outgrid% = 8
    }
    else %outgrid% = 100
/*   }
  end
list %outgrid%.vat
&describe %outgrid%
&type %GRD$MEAN%
&if [show display] cn '9999' &then; gridpa %outgrid%
kill %tempgrid%
&if %GRD$MEAN% = 100  or %GRD$MEAN% = %oldmean%  &then; &goto done
&s oldmean = %GRD$MEAN%
rename %outgrid% %tempgrid%
&s oldmean %GRD$MEAN%
&s its [calc %its% + 1]
&goto bigloop

&label done
&type process complete in %its% iterations. %outgrid% is ready.