Changeset 28062

Timestamp:

01/23/24 13:25:41 (14 months ago)

Author:

Mathieu Morlighem

Message:

NEW: adding vz to SIA

File:

• issm/trunk-jpl/src/m/contrib/morlighem/sia.m (modified) (2 diffs)

issm/trunk-jpl/src/m/contrib/morlighem/sia.m

@@ -1 @@
-function [velx,vely,vel]=sia(md)
+function [vx, vy, vz, vel]=sia(md)
 %SIA - computation of Shallow Ice velocities
 %
@@ -6 +6 @@
 %
 %   Usage:
-%      [velx,vely,vel]=sia(md)
+%      [vx, vy, vz, vel]=sia(md)
 
-if md.mesh.dimension~=2,
-        error('Only 2d meshes are allowed to compute velocity balances');
+disp('Info: assuming no basal sliding and isothermal ice');
+
+n=3;
+
+rhog = (md.materials.rho_ice*md.constants.g);
+
+if md.mesh.dimension==2,
+        %Get slope
+        [sx,sy,s]=slope(md);
+
+        %Average thickness and B over all elements.
+        summer=[1;1;1];
+        hel=md.geometry.thickness(md.mesh.elements)*summer/3;
+        Bel=md.materials.rheology_B(md.mesh.elements)*summer/3;
+
+        Ael=Bel.^(-3);
+
+        vx=-2*rhog^n*s.^(n-1).*sx.*Ael/(n+2).*hel.^(n+1);
+        vy=-2*rhog^n*s.^(n-1).*sy.*Ael/(n+2).*hel.^(n+1);
+        vz = zeros(size(vx));
+else
+
+        [sx,sy,s]=slope(md);
+        sx = average(md, sx, 3);
+        sy = average(md, sy, 3);
+        surf = sqrt(sx.^2 + sy.^2);
+        s  = md.geometry.surface;
+        z  = md.mesh.z;
+        H  = md.geometry.thickness;
+
+        A = md.materials.rheology_B.^-3;
+
+        vx =-2*rhog^n*surf.^(n-1).*sx.*A/(n+1).*(H.^(n+1) - (s-z).^(n+1));
+        vy =-2*rhog^n*surf.^(n-1).*sy.*A/(n+1).*(H.^(n+1) - (s-z).^(n+1));
+
+        chi = 1-(s-z)./H;
+        vz  = -md.smb.mass_balance.*((n+2)*chi + (1-chi).^(n+2) -1)/(n+1)
 end
 
-%Get slope
-[sx,sy,s]=slope(md);
-
-%Average thickness and B over all elements.
-summer=[1;1;1];
-hel=md.geometry.thickness(md.mesh.elements)*summer/3;
-Bel=md.materials.rheology_B(md.mesh.elements)*summer/3;
-
-Ael=Bel.^(-3);
-
-velx=-2*(md.materials.rho_ice*md.constants.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
-vely=-2*(md.materials.rho_ice*md.constants.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
-vel=sqrt(velx.^2+vely.^2);
+vel=sqrt(vx.^2+vy.^2+vz.^2);