function md=huttercompute(md);
%HUTTERCOMPUTE - compute the velocity field of the model using Hutter's model
%
%   Usage:
%      md=huttercompute(md)

%define global variables
iceglobal

%determine if run is parallel
if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end;

%for now, only serial support is in
if cluster,
	error('icediagnostic error message: parallel support not implemented yet');
end

%setup some parameters to be passed to the core solution
params.sparsity=md.sparsity;
params.solver_type=md.solver_type;
params.eps_abs=md.eps_abs;
params.eps_rel=md.eps_rel;
params.debug=md.debug;

%First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model
disp('   processing model');
[elements,grids,loads,constraints,materials,part,tpart]=ModelProcessor(md,'diagnostic_hutter');

%generate degrees of freedom
disp('   generating degrees of freedom');
grids=Dof(grids); 

%add constraints to grids. Create y_g vector of constraint values for the entire g set.
disp('   generating constraint vector');
[grids y_g]=SpcGrids(grids,constraints);

%generate grid sets and reduce y_g to ys from g-set to s-set.
disp('   generating grid sets');
BuildGridSets(grids);gridset=gridset; ys=Reducevector_g(y_g);

%Compute velocity, in 3d
gridset=gridset;
slopesurface=zeros(gridset.gsize,1);
slopesurface(1:6:gridset.gsize,1)=md.surface_slopex;
slopesurface(2:6:gridset.gsize,1)=md.surface_slopey;
disp(sprintf('\n%s',['computing Hutter''s velocities...']));
u_g=icediagnostic_core_linear(elements,grids,materials,loads,params,ys,'diagnostic_hutter',struct('slopesurface',slopesurface));
	
md.vx=u_g(1:6:end)*md.yts;
md.vy=u_g(2:6:end)*md.yts;
md.vel=sqrt(md.vx.^2+md.vy.^2);
