Matlab file

%Paths to datasets
modeldatapath='/u/astrid-r1b/larour/ModelData';
thicknesspath=[modeldatapath '/ASEThickness/ASE_thickness'];
surfacepath=[modeldatapath '/BamberDEMAntarctica1km/surface_smooth30_lowslope'];
firnpath=[modeldatapath '/BroekeFirn1km/firn.mat'];

%Read Surface and interpolate onto model
md.surface=InterpFromFile(md.x,md.y,surfacepath,10);

%Correct surface using firn layerfrom Broeke
md.firn_layer=InterpFromFile(md.x,md.y,firnpath,0);
rho_ice=917;
rho_firn=830;
firn_layer_correction=md.firn_layer*(1-rho_firn/rho_ice);
md.surface=md.surface-firn_layer_correction;

%Load thickness and interpolate onto model
load(thicknesspath);
md.thickness=InterpFromGridToMesh(x_m,y_m,thickness,md.x,md.y,10);

%Correct thickness using Bamber's dem
md=ThicknessCorrection(md,15000);

ThicknessCorrection is the following routine:

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;