function [sx,sy,smag]=basalstress(md)
%BASALSTRESS -  evaluates the basal stress
%
%   The basal stress is computed according to the following formula: 
%   sigma=drag^2*Neff^r*vel^s
%   where:
%    Neff is the effective pressure (Neff=rho_ice*g*thickness+rho_water*g*bed),
%    r=q/p
%    s=1/p
%    p and q being friction coefficients found in the literature (Paterson).
%
%   Usage:
%      [Fx,Fy,Fmag]=basalstress(md)

%Compute effective pressure
Neff=md.rho_ice*md.g*md.thickness+md.rho_water*md.g*md.bed;
Neff(find(Neff<0))=0.1;
r=md.q./md.p;
s=1./md.p;
dragel=(md.drag(md.elements(:,1))+md.drag(md.elements(:,2))+md.drag(md.elements(:,3)))/3;
Neffel=(Neff(md.elements(:,1))+Neff(md.elements(:,2))+Neff(md.elements(:,3)))/3;
velel=(md.vel(md.elements(:,1))+md.vel(md.elements(:,2))+md.vel(md.elements(:,3)))/3/md.yts;
vxel=(md.vx(md.elements(:,1))+md.vx(md.elements(:,2))+md.vx(md.elements(:,3)))/3/md.yts;
vyel=(md.vy(md.elements(:,1))+md.vy(md.elements(:,2))+md.vy(md.elements(:,3)))/3/md.yts;

sx=dragel.^2.*Neffel.^r.*velel.^(s-1).*vxel;
sy=dragel.^2.*Neffel.^r.*velel.^(s-1).*vyel;
smag=sqrt(sx.^2+sy.^2);
