Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/loneedges.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/loneedges.m	(revision 22217)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/loneedges.m	(revision 22217)
@@ -0,0 +1,24 @@
+function edges=loneedges(md)
+
+	edges=[];
+
+	for e=1:md.mesh.numberofelements, 
+		for j=1:3, 
+			if j==3,
+				ind1=md.mesh.elements(e,3);
+				ind2=md.mesh.elements(e,1);
+			else
+				ind1=md.mesh.elements(e,j);
+				ind2=md.mesh.elements(e,j+1);
+			end
+
+			%edge ind1 and ind2: 
+			els1=md.mesh.vertexconnectivity(ind1,1: md.mesh.vertexconnectivity(ind1,end));
+			els2=md.mesh.vertexconnectivity(ind2,1: md.mesh.vertexconnectivity(ind2,end));
+			els=intersect(els1,els2);
+			if length(els)~=2,
+				edges=[edges;[ind1,ind2]];
+			end
+			end
+		end
+	end
Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshintersect.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshintersect.m	(revision 22216)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshintersect.m	(revision 22217)
@@ -22,5 +22,10 @@
 		if length(s)>1,
 			indices(i)=s(1);
-			warn('one or more vertices on the global mesh were duplicated!');
+			format long
+			[lats(i) longs(i)]
+			for j=1:length(s),
+				[lat(s(j)) long(s(j))]
+			end
+			error('one or more vertices on the global mesh were duplicated!');
 		else
 			indices(i)=s;
Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshintersect3d.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshintersect3d.m	(revision 22217)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshintersect3d.m	(revision 22217)
@@ -0,0 +1,33 @@
+function indices=meshintersect3d(x,y,z,xs,ys,zs,varargin)
+%MESHINTERSECT - return indices (into x,y and z) of common values between (x,y,z) and (xs,ys,zs). 
+%  i.e: x(index)=xs; y(index)=ys;
+%
+
+
+	%process options: 
+	options=pairoptions(varargin{:});
+
+	%retrieve tolerance: 
+	tolerance=getfieldvalue(options,'tolerance',10);
+
+	%go through lats,longs and find within tolerance, the index of the corresponding value in lat,long: 
+	indices=zeros(length(xs),1);
+	
+	for i=1:length(xs),
+		distance=sqrt((x-xs(i)).^2+(y-ys(i)).^2+(z-zs(i)).^2);
+		s=find(distance<tolerance);
+		if length(s)>1,
+			for j=1:length(s),
+				hold on;plot3(x(s(j)),y(s(j)),z(s(j)),'c.','MarkerSize',40)
+			end
+			distance(s)
+			error(sprintf('one or more vertices on the global mesh were duplicated (offset %i)',i));
+		elseif isempty(s),
+			plot(distance);
+			min(distance);
+			i
+			error('cannot find concurrent vertics!');
+		else
+			indices(i)=s;
+		end
+	end
Index: sm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge.m	(revision 22216)
+++ 	(revision )
@@ -1,89 +1,0 @@
-function md=modelmerge(md1,md2)
-%MODELMERGE  - merge two models by merging their meshes
-%
-%   Usage:
-%      md=modelmerge(md1,md2);
-
-	if nargin~=2 
-		help meshconvert
-		error('meshconvert error message: bad usage');
-	end
-	
-	md=model();
-
-	%first ,copy md1 mesh into md.mesh to initialize: 
-	md.mesh=md1.mesh;
-
-	%some initializatoin: 
-	elements1=md1.mesh.elements;
-	x1=md1.mesh.x;
-	y1=md1.mesh.y;
-	nods1=md1.mesh.numberofvertices;
-	nel1=md1.mesh.numberofelements;
-
-	elements2=md2.mesh.elements;
-	x2=md2.mesh.x;
-	y2=md2.mesh.y;
-	nods2=md2.mesh.numberofvertices;
-	nel2=md2.mesh.numberofelements;
-	segs2=md2.mesh.segments;
-
-	%offset elements2 by nods1: 
-	elements2=elements2+nods1;
-
-	tolerance=1e-5;
-	%go into the vertices on boundary of mesh 1, and figure out which ones are common to mesh2: 
-	verticesonboundary=find(md1.mesh.vertexonboundary); 
-	for i=1:length(verticesonboundary),
-		node1=verticesonboundary(i); xnode1=x1(node1); ynode1=y1(node1);
-		%is there another node with these coordinates in mesh2? 
-		ind=find(sqrt(((x2-xnode1).^2+(y2-ynode1).^2))<tolerance);
-		if ~isempty(ind),
-			x2(ind)=NaN;
-			y2(ind)=NaN;
-			pos=find(elements2==(ind+nods1)); elements2(pos)=node1;
-		end
-	end
-
-	%go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped: 
-	while( ~isempty(find(isnan(x2)))),
-		for i=1:length(x2),
-			if isnan(x2(i)),
-				pos=find(elements2>(i+nods1));
-				elements2(pos)=elements2(pos)-1;
-				x2(i)=[];
-				y2(i)=[];
-				break;
-			end
-		end
-	end
-
-	%merge elements: 
-	elements=[elements1;elements2];
-
-	%merge vertices: 
-	x=[x1;x2]; 
-	y=[y1;y2];
-
-	%output: 
-	md.mesh.x=x;
-	md.mesh.y=y;
-	md.mesh.elements=elements;
-	md.mesh.numberofvertices=length(x);
-	md.mesh.numberofelements=size(elements,1);
-
-	%connectivities: 
-	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
-	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-
-	%find segments: 
-	md.mesh.segments=findsegments(md);
-
-	%vertex on boundary: 
-	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
-	md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-
-	%some checks: 
-	if max(md.mesh.elements)>md.mesh.numberofvertices, 
-		error('issue in modelmerge, one of the element ids is > number of vertices!');
-	end
Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge2d.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge2d.m	(revision 22217)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge2d.m	(revision 22217)
@@ -0,0 +1,88 @@
+function md=modelmerge2d(md1,md2,varargin)
+%MODELMERGE  - merge two models by merging their meshes
+%
+%   Usage:
+%      md=modelmerge(md1,md2);
+	
+	%process options: 
+	options=pairoptions(varargin{:});
+	
+	tolerance=getfieldvalue(options,'tolerance',10^-5);
+	
+	md=model();
+
+	%first ,copy md1 mesh into md.mesh to initialize: 
+	md.mesh=md1.mesh;
+
+	%some initializatoin: 
+	elements1=md1.mesh.elements;
+	x1=md1.mesh.x;
+	y1=md1.mesh.y;
+	nods1=md1.mesh.numberofvertices;
+	nel1=md1.mesh.numberofelements;
+
+	elements2=md2.mesh.elements;
+	x2=md2.mesh.x;
+	y2=md2.mesh.y;
+	nods2=md2.mesh.numberofvertices;
+	nel2=md2.mesh.numberofelements;
+	segs2=md2.mesh.segments;
+
+	%offset elements2 by nods1: 
+	elements2=elements2+nods1;
+
+	%go into the vertices on boundary of mesh 1, and figure out which ones are common to mesh2: 
+	verticesonboundary=find(md1.mesh.vertexonboundary); 
+	for i=1:length(verticesonboundary),
+		node1=verticesonboundary(i); xnode1=x1(node1); ynode1=y1(node1);
+		%is there another node with these coordinates in mesh2? 
+		ind=find(sqrt(((x2-xnode1).^2+(y2-ynode1).^2))<tolerance);
+		if ~isempty(ind),
+			x2(ind)=NaN;
+			y2(ind)=NaN;
+			pos=find(elements2==(ind+nods1)); elements2(pos)=node1;
+		end
+	end
+
+	%go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped: 
+	while( ~isempty(find(isnan(x2)))),
+		for i=1:length(x2),
+			if isnan(x2(i)),
+				pos=find(elements2>(i+nods1));
+				elements2(pos)=elements2(pos)-1;
+				x2(i)=[];
+				y2(i)=[];
+				break;
+			end
+		end
+	end
+
+	%merge elements: 
+	elements=[elements1;elements2];
+
+	%merge vertices: 
+	x=[x1;x2]; 
+	y=[y1;y2];
+
+	%output: 
+	md.mesh.x=x;
+	md.mesh.y=y;
+	md.mesh.elements=elements;
+	md.mesh.numberofvertices=length(x);
+	md.mesh.numberofelements=size(elements,1);
+
+	%connectivities: 
+	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+
+	%find segments: 
+	md.mesh.segments=findsegments(md);
+
+	%vertex on boundary: 
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+
+	%some checks: 
+	if max(md.mesh.elements)>md.mesh.numberofvertices, 
+		error('issue in modelmerge, one of the element ids is > number of vertices!');
+	end
Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge3d.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge3d.m	(revision 22217)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge3d.m	(revision 22217)
@@ -0,0 +1,96 @@
+function md=modelmerge3d(md1,md2,varargin)
+%MODELMERGE  - merge two models by merging their meshes
+%
+%   Usage:
+%      md=modelmerge(md1,md2);
+	
+	%process options: 
+	options=pairoptions(varargin{:});
+	
+	tolerance=getfieldvalue(options,'tolerance',10^-5);
+	
+	md=model();md.mask=maskpsl();
+
+	%first ,copy md1 mesh into md.mesh to initialize: 
+	md.mesh=md1.mesh;
+
+	%some initializatoin: 
+	elements1=md1.mesh.elements;
+	x1=md1.mesh.x;
+	y1=md1.mesh.y;
+	z1=md1.mesh.z;
+	nods1=md1.mesh.numberofvertices;
+	nel1=md1.mesh.numberofelements;
+
+	elements2=md2.mesh.elements;
+	x2=md2.mesh.x;
+	y2=md2.mesh.y;
+	z2=md2.mesh.z;
+	nods2=md2.mesh.numberofvertices;
+	nel2=md2.mesh.numberofelements;
+
+	%offset elements2 by nods1: 
+	elements2=elements2+nods1;
+
+	%go into the vertices on boundary of mesh 1, and figure out which ones are common to mesh2: 
+	verticesonboundary=find(md1.mesh.vertexonboundary); 
+	for i=1:length(verticesonboundary),
+		node1=verticesonboundary(i); xnode1=x1(node1); ynode1=y1(node1); znode1=z1(node1);
+		%is there another node with these coordinates in mesh2? 
+		ind=find(sqrt(((x2-xnode1).^2+(y2-ynode1).^2)+(z2-znode1).^2)<tolerance);
+		if length(ind)>1,
+			disp('should reduce the tolerance, several vertices picked up!');
+		end
+		if ~isempty(ind),
+			x2(ind)=NaN;
+			y2(ind)=NaN;
+			z2(ind)=NaN;
+			pos=find(elements2==(ind+nods1)); elements2(pos)=node1;
+		end
+	end
+
+	%go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped: 
+	while( ~isempty(find(isnan(x2)))),
+		for i=1:length(x2),
+			if isnan(x2(i)),
+				pos=find(elements2>(i+nods1));
+				elements2(pos)=elements2(pos)-1;
+				x2(i)=[];
+				y2(i)=[];
+				z2(i)=[];
+				break;
+			end
+		end
+	end
+
+	%merge elements: 
+	elements=[elements1;elements2];
+
+	%merge vertices: 
+	x=[x1;x2]; 
+	y=[y1;y2];
+	z=[z1;z2];
+
+	%output: 
+	md.mesh.x=x;
+	md.mesh.y=y;
+	md.mesh.z=z;
+	md.mesh.elements=elements;
+	md.mesh.numberofvertices=length(x);
+	md.mesh.numberofelements=size(elements,1);
+
+	%connectivities: 
+	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+
+	%find segments: 
+	md.mesh.segments=findsegments(md);
+
+	%vertex on boundary: 
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+
+	%some checks: 
+	if max(md.mesh.elements)>md.mesh.numberofvertices, 
+		error('issue in modelmerge, one of the element ids is > number of vertices!');
+	end