0001 function [velx,vely,vel]=balvel(md)
0002
0003
0004
0005
0006
0007
0008
0009
0010 if ~strcmpi(md.type,'2d'),
0011 error('Only 2d meshes are allowed to compute velocity balances');
0012 end
0013
0014
0015 [sx,sy,s]=slope(md);
0016
0017
0018 summer=[1;1;1];
0019 hel=md.thickness(md.elements)*summer/3;
0020 Bel=md.B(md.elements)*summer/3;
0021
0022 Ael=Bel.^(-3);
0023
0024 velx=-2*(md.rho_ice*md.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
0025 vely=-2*(md.rho_ice*md.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
0026 vel=sqrt(velx.^2+vely.^2);