function [velx,vely,vel]=balvel(md)
%BALVEL - computation of balanced velocities
%
%   This routine uses the model of Hutter to compute the velocities
%   of a 2d model using the surface slope
%
%   Usage:
%      [velx,vely,vel]=balvel(md)

if ~strcmpi(md.type,'2d'),
	error('Only 2d meshes are allowed to compute velocity balances');
end

%Get slope
[sx,sy,s]=slope(md);

%Average thickness and B over all elements.
summer=[1;1;1];
hel=md.thickness(md.elements)*summer/3;
Bel=md.B(md.elements)*summer/3;

Ael=Bel.^(-3);

velx=-2*(md.rho_ice*md.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
vely=-2*(md.rho_ice*md.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
vel=sqrt(velx.^2+vely.^2);
