Index: /issm/trunk-jpl/src/m/mesh/patchglobe.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/patchglobe.m	(revision 20110)
+++ /issm/trunk-jpl/src/m/mesh/patchglobe.m	(revision 20111)
@@ -5,9 +5,11 @@
 
 	%recover basic options:
-	hem=getfieldvalue(options,'hem');
 	bandwidth=getfieldvalue(options,'bandwidth',100000);
 
 	%give ourselves a unique temporary directory: 
 	temproot=tempname; mkdir(temproot);
+
+	%align external segments on 2d model: 
+	mh2d.segments=alignsegments(mh2d.segments);
 
 	%figure out domain outline: 
@@ -18,11 +20,7 @@
 	expcontract([temproot '/PatchBroad.exp'],[temproot '/PatchBroad.exp'],-bandwidth);
 
-	%now flag vertices that are within the contour:
-	[x,y]=ll2xy(mh.lat,mh.long,hem);
-	if hem==1,
-		flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1) & mh.lat>=0;
-	else
-		flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1) & mh.lat<=0;
-	end
+	%now flag vertices (from mh2d's broad contour that are on the global mesh: do this in the local 2d mesh reference system. 
+	[x,y]=gdaltransform(mh.long,mh.lat,'EPSG:4326',['EPSG:' num2str(mh2d.epsg)]);
+	flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1);
 
 	%expand flags to any element that touches the contour: 
@@ -37,12 +35,14 @@
 
 	%x,y for segments: 
-	[xsegs,ysegs]=ll2xy(mh.lat(mh.segments(:,1)),mh.long(mh.segments(:,1)),hem);
+	[xsegs,ysegs]=gdaltransform(mh.long(mh.segments(:,1)),mh.lat(mh.segments(:,1)),'EPSG:4326',['EPSG:' num2str(mh2d.epsg)]);
 
-	%create contour out of these segments:
+	%create lat,long contour out of these segments:
 	meshtodomain(mh,[temproot '/PatchEnveloppe.exp'],'latlong','on');
-	expll2xy([temproot '/PatchEnveloppe.exp'],hem);
+		
+	%get these lat,long transformed to local mesh referencial:
+	env=expread([temproot '/PatchEnveloppe.exp']); 
+	[env.x,env.y]=gdaltransform(env.x,env.y,'EPSG:4326',['EPSG:' num2str(mh2d.epsg)]);
 
 	%now, create domain outine from broad enveloppe and initial mesh 
-	env=expread([temproot '/PatchEnveloppe.exp']); 
 	dom=expread([temproot '/Patch.exp']);
 	
@@ -61,4 +61,7 @@
 	expwrite(domain,[temproot '/PatchBand.exp']);
 
+	%plot mesh: 
+	if  getfieldvalue(options,'plot',0), expdisp([temproot '/PatchBand.exp']); end
+
 	%mesh: 
 	mdb=bamg(model(),'domain',[temproot '/PatchBand.exp'],'MaxCornerAngle',1e-15,'gradation',10000); 
@@ -71,5 +74,5 @@
 
 	%augment patch with band
-	[mhb.lat,mhb.long]=xy2ll(mhb.x,mhb.y,hem);
+	[mhb.long,mhb.lat]=gdaltransform(mhb.x,mhb.y,['EPSG:' num2str(mh2d.epsg)],'EPSG:4326');
 	mh2db=augment2dmesh(mh2d,mhb);