Index: /issm/trunk/src/m/classes/public/ThicknessCorrection.m
===================================================================
--- /issm/trunk/src/m/classes/public/ThicknessCorrection.m	(revision 855)
+++ /issm/trunk/src/m/classes/public/ThicknessCorrection.m	(revision 856)
@@ -1,39 +1,60 @@
-function md=ThicknessCorrection(md,filename,distance)
+function md=ThicknessCorrection(md,varargin)
 %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)
+%   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,filename,distance)
+%      md=ThicknessCorrection(md,varargin);
+%
+%   Example:
+%      md=ThicknessCorrection(md);
+%      md=ThicknessCorrection(md,15000);
 
-%some checks
-if ~exist(filename),
-	error(['ThicknessCorrection error message: the file ' filename  ' does not exist']);
+%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),:)));
+
+%get distance
+if nargin==2,
+	distance=varargin{1};
+else
+	distance=10000;
 end
 
-in=ContourToMesh(md.elements,md.x,md.y,expread(filename,1),'node',1);
-pos_shelf=find(in & ~md.gridonicesheet);
-pos_sheet=find(in & md.gridonicesheet);
+%modify thickness
 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)));
+	[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
+		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
-di=md.rho_ice/md.rho_water;
-minth=1/(1-di);
+minth=1/(1-md.rho_ice/md.rho_water);
 pos=find(isnan(md.thickness) | (md.thickness<=0));
-md.thickness(pos)=minth;
+thickness(pos)=minth;
 
 %change bed to take into account the changes in thickness
+md.thickness=thickness;
 md.bed=md.surface-md.thickness;