g.region -d

r.mapcalc --overwrite "dem = 0.5*x()" #elevation map
r.mapcalc --overwrite "lithology = 1" #lithology map (defining the soil class for those areas where no layers are defined)
r.mapcalc --overwrite "depth1 = max(0,-10+x()/10)" #depth of layer bottom maps
r.mapcalc --overwrite "depth2 = max(0,-5+x()/10)"
r.mapcalc --overwrite "depth3 = max(0,0+x()/10)"
r.mapcalc --overwrite "depth4 = max(0,5+x()/10)"
r.mapcalc --overwrite "depth5 = max(0,20+x()/10)"
r.mapcalc --overwrite "gwdepth = max(0,0.5*x()-40)" #depth of groundwater table

#Soil class-based mode (one soil class with five layers)
r.slope.stability -t -v prefix=so_sing_class cellsize=1 model=c elevation=dem soilclass=lithology gwdepth=gwdepth seepage=1 numlayers=5 depthmaps=depth1,depth2,depth3,depth4,depth5 geotech=1,1,15000,4000,30,40,1,2,15000,2000,20,35,1,3,15000,2500,40,40,1,4,15000,5000,30,40,1,5,27000,100000,89,0 elldens=0 ellips=100,100,50,50,25,25,0,0,3,1,100,100,-9999,-9999

#Layer-based mode (five layers, each associated to one of four soil classes)
r.slope.stability -t -v prefix=so_sing_lyr cellsize=1 model=lr elevation=dem soilclass=lithology gwdepth=gwdepth seepage=2 prelayers=depth classlayers=2,1,3,1,4 geotech=1,15000,4000,20,40,2,15000,5000,25,35,3,15000,2500,40,40,4,27000,100000,89,0 elldens=0 ellips=100,100,50,50,25,25,0,0,3,1,100,100,-9999,-9999 maxlayers=6

g.region -d