Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 19954)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 19955)
@@ -303,5 +303,5 @@
 			md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
 			md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-			md.mesh.segments=contourenvelope(md);
+			md.mesh.segments=contourenvelope(md.mesh);
 
 		end % }}}
Index: /issm/trunk-jpl/src/m/exp/meshtodomain.m
===================================================================
--- /issm/trunk-jpl/src/m/exp/meshtodomain.m	(revision 19955)
+++ /issm/trunk-jpl/src/m/exp/meshtodomain.m	(revision 19955)
@@ -0,0 +1,39 @@
+function meshtodomain(mh,domainname,varargin)
+%MESHTODOMAIN - recover a domain outline  from a model's mesh's segments
+%
+%   Usage:
+%      meshtodomain(mh,domainname,varargin)
+%
+%   Example:
+%      meshtodomain(md.mesh,'domainoutline.exp','latlong','on');
+%
+%   See also EXPREAD
+
+	%handle options: 
+	options=pairoptions(varargin{:});
+
+	segments=mh.segments; nt=length(segments);
+
+	%build domain contour: 
+	x=[];  y=[]; 
+
+	if strcmpi(getfieldvalue(options,'latlong','off'),'on'),
+		if isnan(mh.lat) | isnan(mh.long), error('meshtodomain error message: requested domain be output in lat,long referential, but mesh does not contain this information!'); end
+		for i=1:nt,
+		   x=[x;mh.long(segments(i,1))];
+		   y=[y;mh.lat(segments(i,1))];
+		end
+	else
+		for i=1:nt,
+		   x=[x;mh.x(segments(i,1))];
+		   y=[y;mh.y(segments(i,1))];
+		end
+	end
+
+	domain.x=x; 
+	domain.y=y; 
+	domain.density=1; 
+	domain.name=domainname;
+
+	expwrite(domain,domainname);
+end
Index: /issm/trunk-jpl/src/m/miscellaneous/alignsegments.m
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/alignsegments.m	(revision 19955)
+++ /issm/trunk-jpl/src/m/miscellaneous/alignsegments.m	(revision 19955)
@@ -0,0 +1,19 @@
+function newsegments=alignsegments(segments)
+%ALIGNSEGMENTS: 
+% 
+%
+
+	nt=length(segments);
+	newsegments=zeros(nt,3);
+	newsegments(1,:)=segments(1,:);
+
+	for  i=2:nt, 
+		last=newsegments(i-1,2); %last vertex of the previous segment: 
+		for j=1:nt,
+			if last==segments(j,1),
+				%we found the next segment: 
+				newsegments(i,:)=segments(j,:);
+				break;
+			end
+		end
+	end
Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 19954)
+++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 19955)
@@ -1,11 +1,11 @@
-function segments=contourenvelope(md,varargin)
+function segments=contourenvelope(mh,varargin)
 %CONTOURENVELOPE - build a set of segments enveloping a contour .exp
 %
 %   Usage:
-%      segments=contourenvelope(md,varargin)
+%      segments=contourenvelope(mh,varargin)
 %
 %   Example:
-%      segments=contourenvelope(md,'Stream.exp');
-%      segments=contourenvelope(md);
+%      segments=contourenvelope(mh,'Stream.exp');
+%      segments=contourenvelope(mh);
 
 %some checks
@@ -33,25 +33,25 @@
 %Now, build the connectivity tables for this mesh.
 %Computing connectivity
-if (size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices & size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices2d),
-	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+if isnan(mh.vertexconnectivity),
+	mh.vertexconnectivity=NodeConnectivity(mh.elements,mh.numberofvertices);
 end
-if (size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements & size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements2d),
-	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+if isnan(mh.elementconnectivity),
+	mh.elementconnectivity=ElementConnectivity(mh.elements,mh.vertexconnectivity);
 end
 
 %get nodes inside profile
-mesh.elementconnectivity=md.mesh.elementconnectivity;
-if dimension(md.mesh)==2,
-	mesh.elements=md.mesh.elements;
-	mesh.x=md.mesh.x;
-	mesh.y=md.mesh.y;
-	mesh.numberofvertices=md.mesh.numberofvertices;
-	mesh.numberofelements=md.mesh.numberofelements;
+mesh.elementconnectivity=mh.elementconnectivity;
+if dimension(mh)==2,
+	mesh.elements=mh.elements;
+	mesh.x=mh.x;
+	mesh.y=mh.y;
+	mesh.numberofvertices=mh.numberofvertices;
+	mesh.numberofelements=mh.numberofelements;
 else
-	mesh.elements=md.mesh.elements2d;
-	mesh.x=md.mesh.x2d;
-	mesh.y=md.mesh.y2d;
-	mesh.numberofvertices=md.mesh.numberofvertices2d;
-	mesh.numberofelements=md.mesh.numberofelements2d;
+	mesh.elements=mh.elements2d;
+	mesh.x=mh.x2d;
+	mesh.y=mh.y2d;
+	mesh.numberofvertices=mh.numberofvertices2d;
+	mesh.numberofelements=mh.numberofelements2d;
 end