function dhdt=thicknessevolution(md)
%THICKNESSEVOLUTION - compute the new thickness of a model after ∆t
%
%   This routine compute the new thickness of a model after a time step
%   according to the following formula:
%   dh/dt=div(Hu)+Ms-Mb
%
%   Usage:
%      dhdt=thicknessevolution(md)

if (length(md.vx)~=md.numberofgrids)|(length(md.vy)~=md.numberofgrids)
	error('thicknessevolution error message: velocity empty')
end

alpha=zeros(md.numberofelements,3);
beta=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,:);
end

summation=[1;1;1];
dHu=(md.vx(md.elements).*md.thickness(md.elements).*alpha)*summation;
dHv=(md.vy(md.elements).*md.thickness(md.elements).*beta)*summation;

dHu=averaging(md,dHu,0);
dHv=averaging(md,dHv,0);

dhdt=dHu+dHv+md.accumulation-md.melting;
