function md=ThicknessCorrection(md,filename,distance)
%THICKNESSCORRECTION - correct the thickness of the ice shelf near the grounding line
%
%   This routine corrects the thickness and the bed at the transition ice sheet - shelf
%   due to the assumption of hydrostatic equilibrium.
%   it treats the area given by the Argus file given in input as follows:
%   1: linearize the transition iceshelf ice sheet
%    thickness = coeff * thickness_GL + (1-coeff) * thickness_shelf
%    coeff = min(0,ditance to GL / distance)
%   2: take the minimum between this linearized thickness and the previous thickness
%    thickness = min(linearized thickness, previous thickness)
%
%   Usage:
%      md=ThicknessCorrection(md,filename,distance)

%some checks
if ~exist(filename),
	error(['ThicknessCorrection error message: the file ' filename  ' does not exist']);
end

in=ArgusContourToMesh(md.elements,md.x,md.y,expread(filename,1),'node',1);
pos_shelf=find(in & md.gridoniceshelf);
pos_sheet=find(in & ~md.gridoniceshelf);
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_sheet)).^2+(md.y(pos_shelf(i))-md.y(pos_sheet)).^2));
	%thickness is minimum between hydrostatic equilibrium and thickness of the closest grid on ice sheet
	coeff=min(1,d/distance);
	md.thickness(pos_shelf(i))=min(coeff*md.thickness(pos_shelf(i))+(1-coeff)*md.thickness(pos_sheet(posd)),md.thickness(pos_shelf(i)));
end

%check the computed thickness
di=md.rho_ice/md.rho_water;
minth=1/(1-di);
pos=find(isnan(md.thickness) | (md.thickness<=0));
md.thickness(pos)=minth;

%change bed to take into account the changes in thickness
md.bed=md.surface-md.thickness;
