Index: /issm/trunk-jpl/src/m/parameterization/killberg.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/killberg.m	(revision 24110)
+++ /issm/trunk-jpl/src/m/parameterization/killberg.m	(revision 24110)
@@ -0,0 +1,60 @@
+%function ice_levelset = killberg(md)
+function ice_levelset = killberg(md)
+%KILLBERG - kill ice berg
+%
+%   check md.mask.ice_levelset and md.mask.groundedice_levelset and
+%   remove icebergs from md.mask.ice_levelset. This function is adapted
+%   from src/c/modules/KillIcebergsx/KillIcebergsx.cpp
+%
+%   Usage:
+%      ice_levelset = killberg(md)
+
+%Get number of vertices per element
+nbv_per_element = size(md.mesh.elements,2);
+
+%Initialize vectors (mask = is active and connected to grounded ice)
+mask         = zeros(md.mesh.numberofvertices,1);
+element_flag = zeros(md.mesh.numberofelements,1);
+
+%Step 1, go through all elements and put 1 in local_mask if the element is grounded
+disp('Looking for isolated patches of floating ice (icebergs)');
+isice      = min(md.mask.ice_levelset(md.mesh.elements),[],2)<0;
+isgrounded = max(md.mask.groundedice_levelset(md.mesh.elements),[],2)>0;
+pos = find(isgrounded);
+element_flag(pos) = 1;
+mask(md.mesh.elements(pos,:)) = 1;
+
+iter = 1;
+more = true;
+while(more)
+	disp(['   -- iteration ' num2str(iter)]);
+	more = false;
+
+	pos = find(~element_flag & isice);
+	if numel(pos)==0
+		break;
+	end
+
+	for i=pos'
+		indices = md.mesh.elements(i,:);
+		MAX=max(mask(indices));
+		if(MAX==0)
+			continue;
+		else
+			element_flag(i) = 1;
+			mask(indices) = 1;
+			more = true;
+		end
+	end
+	iter = iter+1;
+end
+
+%OK now change ice_levelset accroding to mask
+pos = find(mask==0 & md.mask.ice_levelset<0);
+if numel(pos)
+	disp(['REMOVING ' num2str(numel(pos)) ' vertices on icebergs']);
+	ice_levelset = md.mask.ice_levelset;
+	ice_levelset(pos) = +1;
+else
+	disp('No iceberg found!');
+end