md=triangle(model(),'../Exp/Square.exp',150000.);
md=setmask(md,'','');
md=parameterize(md,'../Par/SquareSheetConstrained.par');
md=setflowequation(md,'SSA','all');
md.cluster=generic('name',oshostname(),'np',3);
md.initialization.watercolumn=10^-3*ones(md.mesh.numberofvertices,1);
md.hydrology.spcwatercolumn=NaN*ones(md.mesh.numberofvertices,1);
pos=find(md.mesh.y==0);
md.hydrology.spcwatercolumn(pos)=2*10^-3;
md=solve(md,HydrologySolutionEnum());

%Fields and tolerances to track changes
field_names     ={'HydrologyWaterVx1','HydrologyWaterVy1','WaterColumn1', 'HydrologyWaterVx2','HydrologyWaterVy2','WaterColumn2', 'HydrologyWaterVx3','HydrologyWaterVy3','WaterColumn3'};
field_tolerances={1e-09,1e-01,1e-09, 1e-09,1e-01,1e-09, 1e-09,1e-01,1e-09};
field_values={...
	(md.results.HydrologySolution(1).HydrologyWaterVx),...
	(md.results.HydrologySolution(1).HydrologyWaterVy),...
	(md.results.HydrologySolution(1).Watercolumn),...
	(md.results.HydrologySolution(2).HydrologyWaterVx),...
	(md.results.HydrologySolution(2).HydrologyWaterVy),...
	(md.results.HydrologySolution(2).Watercolumn),...
	(md.results.HydrologySolution(3).HydrologyWaterVx),...
	(md.results.HydrologySolution(3).HydrologyWaterVy),...
	(md.results.HydrologySolution(3).Watercolumn),...
	};
