Changeset 1314

Timestamp:

07/14/09 14:57:42 (16 years ago)

Author:

Eric.Larour

Message:

improved basal stress computation

File:

• issm/trunk/src/m/classes/public/basalstress.m (modified) (1 diff)

issm/trunk/src/m/classes/public/basalstress.m

@@ -1 @@
-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).
+function [bx by b]=basalstress(md,options_structure,width,i,type);
+%BASALSTRESS - compute basal stress from basal drag and geometric information. 
 %
 %   Usage:
-%      [Fx,Fy,Fmag]=basalstress(md)
+%      basalstress(md,options_structure,width,i,type);
+%
+%   See also: plot_basaldrag
 
-%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;
+%check layer
+if strcmpi(md.type,'3d')
+        if options_structure.layer~=1,
+                disp('plot_basaldrag warning: basal drag is displayed in the lower layer')
+                options_structure.layer=1;
+        end
+end
 
-sx=dragel.^2.*Neffel.^r.*velel.^(s-1).*vxel;
-sy=dragel.^2.*Neffel.^r.*velel.^(s-1).*vyel;
-smag=sqrt(sx.^2+sy.^2);
+%compute exponents
+s=averaging(md,1./md.p,0);
+r=averaging(md,md.p./md.q,0);
+
+%compute horizontal velocity
+ub=sqrt(md.vx.^2+md.vy.^2)/md.yts;
+ubx=md.vx/md.yts;
+uby=md.vy/md.yts;
+
+%compute basal drag
+bx=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag).^2.*ubx.^s;
+by=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag).^2.*uby.^s;
+b=(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)).^r.*(md.drag).^2.*ub.^s;