function wb=basevert(md)
%BASEVERT - computes the basal vertical velcities
%
%   This routine computes the basal vertical velocities of ice shelves
%   for 2d models only using the following formula:
%   wb=rho_ice/rho_water*div(thickness*vel_horiz)+vel_horiz.grad(base)
%
%   Usage:
%      wb=basevert(md);

alpha=zeros(md.numberofelements,3);
beta=zeros(md.numberofelements,3);
gamma=zeros(md.numberofelements,3);

for n=1:md.numberofelements
	X=inv([md.x(md.elements(n,:)) md.y(md.elements(n,:)) ones(3,1)]);
	alpha(n,:)=X(1,:);
	beta(n,:)=X(2,:);
	gamma(n,:)=X(3,:);
end

hu=md.thickness.*md.vx;
hv=md.thickness.*md.vy;

summation=[1;1;1];
hux=(hu(md.elements).*alpha)*summation;
hvy=(hv(md.elements).*beta)*summation;

uelem=md.vx(md.elements)*summation/3;
velem=md.vy(md.elements)*summation/3;

dbdx=(md.bed(md.elements).*alpha)*summation;
dbdy=(md.bed(md.elements).*beta)*summation;

wb=-md.rho_ice/md.rho_water*(hux+hvy)+uelem.*dbdx+velem.*dbdy;
