Index: /issm/trunk-jpl/src/m/classes/basin.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/basin.m	(revision 24155)
+++ /issm/trunk-jpl/src/m/classes/basin.m	(revision 24156)
@@ -9,4 +9,5 @@
 		epsg              = 3426;
 		name              = '';
+		continent         = '';
 	end
 	methods (Static)
@@ -29,4 +30,5 @@
 					self.boundaries=getfieldvalue(options,'boundaries',{});
 					self.name=getfieldvalue(options,'name','');
+					self.continent=getfieldvalue(options,'continent','');
 					self.epsg=getfieldvalue(options,'epsg',3426);
 			end
@@ -34,4 +36,5 @@
 		function self = setdefaultparameters(self) % {{{
 			self.name='';
+			self.continent='';
 			self.epsg=3426;
 			self.boundaries={};
@@ -40,4 +43,5 @@
 		function disp(self) % {{{
 			disp(sprintf('   basin parameters:'));
+			fielddisplay(self,'continent','continent name');
 			fielddisplay(self,'name','basin name');
 			fielddisplay(self,'epsg','epsg projection number for the entire basin');
@@ -47,4 +51,22 @@
 			end
 
+		end % }}}
+		function boolean=isnameany(self,varargin) % {{{
+			boolean=0;
+			for  i=1:length(varargin),
+				if strcmpi(self.name,varargin{i}), 
+					boolean=1;
+					break;
+				end
+			end
+		end % }}}
+		function boolean=iscontinentany(self,varargin) % {{{
+			boolean=0;
+			for  i=1:length(varargin),
+				if strcmpi(self.continent,varargin{i}), 
+					boolean=1;
+					break;
+				end
+			end
 		end % }}}
 		function output=outputname(self,varargin) % {{{
@@ -87,10 +109,16 @@
 				boundary=self.boundaries{i};
 				contour=boundary.edges();
-				[contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.projection()),sprintf('EPSG:%i',self.epsg));
+				[contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.epsg),sprintf('EPSG:%i',self.epsg));
 				x=[x;contour.x];
 				y=[y;contour.y];
 			end
+			%close the contour: 
+			if x(end)~=x(1) | y(end)~=y(1), 
+				x(end)=x(1); y(end)=y(1);
+			end
+
 			out.x=x;
 			out.y=y;
+			out.nods=length(x);
 		end % }}}
 		function output=shapefilecrop(self,varargin) % {{{
Index: /issm/trunk-jpl/src/m/classes/boundary.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/boundary.m	(revision 24155)
+++ /issm/trunk-jpl/src/m/classes/boundary.m	(revision 24156)
@@ -66,7 +66,4 @@
 		end
 		end % }}}
-		function output=projection(self) % {{{
-			output=self.epsg;
-		end % }}}
 		function plot(self,varargin) % {{{
 			%recover options
Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh3dsurface.m	(revision 24155)
+++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.m	(revision 24156)
@@ -17,4 +17,5 @@
 		long                        = NaN;
 		r                           = NaN;
+		area                        = NaN;
 
 		vertexonboundary            = NaN;
@@ -115,4 +116,5 @@
 			fielddisplay(obj,'long','vertices longitude [degrees]');
 			fielddisplay(obj,'r','vertices radius [m]');
+			fielddisplay(obj,'area','elemental areas [m^2]');
 			
 			fielddisplay(obj,'edges','edges of the 2d mesh (vertex1 vertex2 element1 element2)');
@@ -180,4 +182,5 @@
 			writejs1Darray(fid,[modelname '.mesh.long'],self.long);
 			writejs1Darray(fid,[modelname '.mesh.r'],self.r);
+			writejs1Darray(fid,[modelname '.mesh.area'],self.area);
 			writejs1Darray(fid,[modelname '.mesh.vertexonboundary'],self.vertexonboundary);
 			writejs2Darray(fid,[modelname '.mesh.edges'],self.edges);
Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 24155)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 24156)
@@ -9,9 +9,10 @@
 %      slm = sealevel('icecap',md_greenland,'icecap',md_antarctica,'earth',md_earth);
 
-classdef sealevelmodel
+classdef sealevelmodel < handle
 	properties (SetAccess=public) %Model fields
 		% {{{
 		icecaps          = {}; % list of ice cap models
 		earth            = 0;  % model for the whole earth
+		basins           = {}; % list  of basins, matching icecaps, where shapefile info is held.
 		cluster          = 0;
 		miscellaneous    = 0;
@@ -19,5 +20,6 @@
 		private          = 0;
 		range            = 0;
-		mergedcaps     = 0;
+		mergedcaps       = 0;
+		transitions      = {};
 		%}}}
 	end
@@ -87,4 +89,5 @@
 			slm.cluster           = generic();
 			slm.range             = {};
+			slm.transitions       = {};
 		end
 		%}}}
@@ -128,5 +131,5 @@
 							continue;
 						end
-											end
+					end
 				end
 				self.mergedcaps{2*(i-1)+1}=md;
@@ -137,4 +140,177 @@
 				disp(sprintf('%i: %s',i,self.icecaps{i}.miscellaneous.name));
 			end
+		end % }}}
+		function addbasin(self,bas) % {{{
+		if ~strcmpi(class(bas),'basin')
+			error('addbasin method only takes a ''basin'' class object as input');
+		end;
+		self.basins{end+1}=bas;
+		end % }}}
+		function intersections(self,varargin) % {{{
+
+		options=pairoptions(varargin{:}); 
+		force=getfieldvalue(options,'force',0);
+
+		for i=1:length(self.icecaps),
+			mdi=self.icecaps{i};
+			mdi=TwoDToThreeD(mdi);
+		
+			disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+		
+			self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
+
+		end
+
+		end % }}}
+		function checkintersections(self) % {{{
+		flags=zeros(self.earth.mesh.numberofvertices,1);
+		for i=1:length(self.basins),
+			flags(self.transitions{i})=i;
+		end
+		plotmodel(self.earth,'data',flags,'coastline','on');
+
+		end % }}}
+		function baslist=basinindx(self,varargin) % {{{
+			options=pairoptions(varargin{:}); 
+			continent=getfieldvalue(options,'continent','all');
+			bas=getfieldvalue(options,'basin','all');
+
+			%expand continent list: {{{
+			if iscell(continent),
+				if length(continent)==1,
+					 if strcmpi(continent{1},'all'),
+						 %need to transform this into a list of continents: 
+						 continent={}; 
+						 for i=1:length(self.basins),
+							 continent{end+1}=self.basins{i}.continent;
+						 end
+						 continent=unique(continent);
+					 end
+				else
+					%nothing to do, we have a list of continents
+				end
+			else
+				if strcmpi(continent,'all'),
+					%need to transform this into a list of continents: 
+					continent={}; 
+					for i=1:length(self.basins),
+						 continent{end+1}=self.basins{i}.continent;
+					end
+					continent=unique(continent);
+				else
+					continent={continent};
+				end
+			end
+			%}}}
+			%expand basins list using the continent list above and the extra bas discriminator: %{{{
+			if iscell(bas),
+				if length(bas)==1,
+					 if strcmpi(bas{1},'all'),
+						 %need to transform this into a list of basins: 
+						 baslist=[];
+						 for i=1:length(self.basins),
+							 if self.basins{i}.iscontinentany(continent{:}),
+								 baslist(end+1)=i;
+							 end
+						 end
+						 baslist=unique(baslist);
+					 else
+						 bas=bas{1};
+						 baslist=[];
+						 for i=1:length(self.basins),
+							 if self.basins{i}.iscontinentany(continent{:}),
+								 if self.basins{i}.isnameany(bas),
+									 baslist(end+1)=i;
+								 end
+							 end
+						 end
+
+					 end
+				else
+					%we have a list of basin names: 
+					baslist=[];
+					for i=1:length(bas),
+						basname=bas{i};
+						for j=1:length(self.basins),
+							if self.basins{j}.iscontinentany(continent{:}),
+								if self.basins{j}.isnameany(basname),
+									baslist(end+1)=j;
+								end
+							end
+						end
+						baslist=unique(baslist);
+					end
+				end
+			else
+				if strcmpi(bas,'all'),
+					baslist=[];
+					for i=1:length(self.basins),
+						if self.basins{i}.iscontinentany(continent{:}),
+							baslist(end+1)=i;
+						end
+					end
+					baslist=unique(baslist);
+				else
+					baslist=[];
+					for i=1:length(self.basins),
+						if self.basins{i}.iscontinentany(continent{:}),
+							if self.basins{i}.isnameany(bas),
+								baslist(end+1)=i;
+							end
+						end
+					end
+					baslist=unique(baslist);
+				end
+			end
+			%}}}
+
+		end % }}}
+		function addicecap(self,md) % {{{
+		if ~strcmpi(class(md),'model')
+			error('addicecap method only takes a ''model'' class object as input');
+		end
+		self.icecaps{end+1}=md;
+		end % }}}
+		function basinsplot3d(self,varargin) % {{{
+		for i=1:length(self.basins),
+			self.basins{i}.plot3d(varargin{:});
+		end
+		end % }}}
+		function caticecaps(self,varargin) % {{{
+			
+			%recover options: 
+			options=pairoptions(varargin{:}); 
+			tolerance=getfieldvalue(options,'tolerance',.65);
+			loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
+	
+			%make 3D model:
+			models=self.icecaps;
+			for i=1:length(models),
+				models{i}=TwoDToThreeD(models{i});
+			end
+			
+			%Plug all models together:
+			md=models{1}; 
+			for i=2:length(models),
+				md=modelmerge3d(md,models{i},'tolerance',tolerance);
+				md.private.bamg.landmask=[md.private.bamg.landmask;models{i}.private.bamg.landmask];
+			end
+
+			%Look for lone edges if asked for it: {{{
+			if loneedgesdetect,
+				edges=loneedges(md);
+				plotmodel(md,'data',md.mask.land_levelset);
+				hold on;
+				for i=1:length(edges),
+					ind1=edges(i,1);
+					ind2=edges(i,2);
+					%plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'r*-');
+					plot3([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],[md.mesh.z(ind1),md.mesh.z(ind2)],'g*-');
+				end
+			end %}}}
+	
+			%Plug into earth: 
+			self.earth=md;
+
 		end % }}}
 		function viscousiterations(self) % {{{