Changes between Version 1 and Version 2 of shelf_thickness

Timestamp:

10/12/10 14:39:19 (14 years ago)

Author:

Mathieu Morlighem

Comment:

Added ThicknessCorrection

shelf_thickness

@@ -1 +1 @@
 == Matlab file ==
+
+{{{
+#!m
 
 %Paths to datasets
@@ -23 +26 @@
 %Correct thickness using Bamber's dem
 md=ThicknessCorrection(md,15000);
+}}}
+ThicknessCorrection is the following routine:
+
+{{{
+#!m
+
+function md=ThicknessCorrection(md,distance)
+%THICKNESSCORRECTION - correct the thickness of the ice shelf near the grounding line
+%
+%   This routine corrects the thickness and the bed on the transition zone
+%   by forcing the hydrostatic equilibrium.
+%   the thickness is modified as follows:
+%      thickness = coeff * thickness_observation + (1-coeff) * thickness_hydrostatic
+%   where:
+%      coeff=(d/distance)^2;
+%      distance=10km by default but can be specified
+%
+%   Usage:
+%      md=ThicknessCorrection(md,distance);
+%
+%   Example:
+%      md=ThicknessCorrection(md,15000);
+
+%initialize thickness with the observations, and get hydrostatic thickness from the dem
+thickness=md.thickness;
+thickness_hydro=md.surface/(1-md.rho_ice/md.rho_water);
+
+%get nodes on ice sheet and on ice shelf
+pos_shelf=find(~md.gridonicesheet);
+pos_GL=intersect(unique(md.elements(find(md.elementonicesheet),:)),unique(md.elements(find(md.elementoniceshelf),:)));
+debug=(length(pos_shelf)>50000);
+
+%modify thickness
+if (debug), fprintf('%s','      correction progress:   0.00 %'); end
+for i=1:length(pos_shelf)
+
+   %search the grid on ice sheet the closest to i
+   [d posd]=min(sqrt((md.x(pos_shelf(i))-md.x(pos_GL)).^2+(md.y(pos_shelf(i))-md.y(pos_GL)).^2));
+
+   if d>distance,
+
+      %if d > 15km, hydrostatique equilibrium
+      thickness(pos_shelf(i))=thickness_hydro(pos_shelf(i));
+
+   else
+
+      %else: quadratic combination of hydrostatic equilibrium and observations (this line might be commented if the original thickness is good close to the GL)
+      coeff=(d/distance)^2;
+      thickness(pos_shelf(i))=(1-coeff)*thickness(pos_shelf(i))+coeff*thickness_hydro(pos_shelf(i));
+
+   end
+end
+
+%check the computed thickness
+minth=1/(1-md.rho_ice/md.rho_water);
+pos=find(isnan(md.thickness) | (md.thickness<=0));
+thickness(pos)=minth;
+
+%change bed to take into account the changes in thickness
+md.thickness=thickness;
+md.bed=md.surface-md.thickness;
+}}}