0001 function [sx,sy,smag]=basalstress(md)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 Neff=md.rho_ice*md.g*md.thickness+md.rho_water*md.g*md.bed;
0017 Neff(find(Neff<0))=0.1;
0018 r=md.q./md.p;
0019 s=1./md.p;
0020 dragel=(md.drag(md.elements(:,1))+md.drag(md.elements(:,2))+md.drag(md.elements(:,3)))/3;
0021 Neffel=(Neff(md.elements(:,1))+Neff(md.elements(:,2))+Neff(md.elements(:,3)))/3;
0022 velel=(md.vel(md.elements(:,1))+md.vel(md.elements(:,2))+md.vel(md.elements(:,3)))/3/md.yts;
0023 vxel=(md.vx(md.elements(:,1))+md.vx(md.elements(:,2))+md.vx(md.elements(:,3)))/3/md.yts;
0024 vyel=(md.vy(md.elements(:,1))+md.vy(md.elements(:,2))+md.vy(md.elements(:,3)))/3/md.yts;
0025
0026 sx=dragel.^2.*Neffel.^r.*velel.^(s-1).*vxel;
0027 sy=dragel.^2.*Neffel.^r.*velel.^(s-1).*vyel;
0028 smag=sqrt(sx.^2+sy.^2);