Index: /issm/trunk-jpl/src/m/classes/basin.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/basin.m	(revision 24778)
+++ /issm/trunk-jpl/src/m/classes/basin.m	(revision 24779)
@@ -7,7 +7,7 @@
 	properties (SetAccess=public) 
 		boundaries        = {};
-		epsg              = 3426;
 		name              = '';
 		continent         = '';
+		proj              = epsg2proj(4326);
 	end
 	methods (Static)
@@ -31,5 +31,22 @@
 					self.name=getfieldvalue(options,'name','');
 					self.continent=getfieldvalue(options,'continent','');
-					self.epsg=getfieldvalue(options,'epsg',3426);
+
+					%figure out projection string: 
+					if exist(options,'epsg'),
+						if exist(options,'proj'),
+							error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
+						end
+						epsg=getfieldvalue(options,'epsg');
+						%convert into proj4: 
+						proj=epsg2proj(epsg);
+					elseif exist(options,'proj'),
+						if exist(options,'epsg'),
+							error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
+						end
+						proj=getfieldvalue(options,'proj');
+					else
+						proj= '+proj=longlat +datum=WGS84 +no_defs';
+					end
+					self.proj=proj;
 			end
 		end % }}}
@@ -37,5 +54,5 @@
 			self.name='';
 			self.continent='';
-			self.epsg=3426;
+			self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
 			self.boundaries={};
 
@@ -45,5 +62,5 @@
 			fielddisplay(self,'continent','continent name');
 			fielddisplay(self,'name','basin name');
-			fielddisplay(self,'epsg','epsg projection number for the entire basin');
+			fielddisplay(self,'proj','basin projection string (follows PROJ.4 standard)');
 			fielddisplay(self,'boundaries','list of boundary objects');
 			for i=1:length(self.boundaries),
@@ -87,5 +104,5 @@
 			%add option: 
 			for i=1:length(self.boundaries),
-				self.boundaries{i}.plot('epsg',self.epsg,varargin{:});
+				self.boundaries{i}.plot('proj',self.proj,varargin{:});
 			end
 
@@ -109,5 +126,5 @@
 				boundary=self.boundaries{i};
 				contour=boundary.edges();
-				[contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.epsg),sprintf('EPSG:%i',self.epsg));
+				[contour.x,contour.y]=gdaltransform(contour.x,contour.y,boundary.proj,self.proj);
 				x=[x;contour.x];
 				y=[y;contour.y];
@@ -122,4 +139,26 @@
 			out.nods=length(x);
 		end % }}}
+		function checkconsistency(self,varargin) % {{{
+		
+			%recover options
+			options=pairoptions(varargin{:});
+
+			%figure out if two boundaries are identical: 
+			for i=1:length(self.boundaries),
+				namei=self.boundaries{i}.shpfilename; 
+				for j=i+1:length(self.boundaries),
+					namej=self.boundaries{j}.shpfilename; 
+					if strcmpi(namei,namej),
+						error(['Basin ' self.name ' has two identical boundaries named ' namei ]);
+					end
+				end
+			end
+
+			%go through boundaries and check their consistency: 
+			for i=1:length(self.boundaries),
+				boundary=self.boundaries{i};
+				boundary.checkconsistency();
+			end
+		end % }}}
 		function contourplot(self,varargin) % {{{
 			contour=self.contour();
@@ -134,5 +173,16 @@
 			inshapefile=getfieldvalue(options,'shapefile');
 			outputshapefile=getfieldvalue(options,'outputshapefile','');
-			epsgshapefile=getfieldvalue(options,'epsgshapefile');
+
+			if (exist(options,'epsgshapefile') &  exist(options,'projshapefile')), 
+				error('Basin shapefilecrop error message: cannot specify epsgshapefile and projshapefile at the same time');
+			end
+			if exist(options,'epsgshapefile'),
+				projshapefile=epsg2proj(getfieldvalue(options,'epsgshapefile'));
+			elseif exist(options,'projshapefile'),
+				projshapefile=getfieldvalue(options,'projshapefile');
+			else
+				error('Basin shapefilecrop error message: epsgshapefile or projshapefile should be specified');
+			end
+
 
 			%create list of contours that have critical length > threshold:  (in lat,long)
@@ -150,11 +200,9 @@
 
 			%project onto reference frame:
-			if self.epsg~=epsgshapefile,
-				for i=1:length(contours),
-					h=contours(i);
-					[h.x,h.y]=gdaltransform(h.x,h.y,sprintf('EPSG:%i',epsgshapefile),sprintf('EPSG:%i',self.epsg));
-					contours(i).x=h.x;
-					contours(i).y=h.y;
-				end
+			for i=1:length(contours),
+				h=contours(i);
+				[h.x,h.y]=gdaltransform(h.x,h.y,projshapefile,self.proj);
+				contours(i).x=h.x;
+				contours(i).y=h.y;
 			end
 
Index: /issm/trunk-jpl/src/m/classes/boundary.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/boundary.m	(revision 24778)
+++ /issm/trunk-jpl/src/m/classes/boundary.m	(revision 24779)
@@ -9,5 +9,6 @@
 		shpfilename       = '';
 		orientation       = 'normal';  %other value is 'reverse'
-		epsg              = 4326; %EPSG number, default value is WGS 84 Lat,Long reference frame.
+		proj              = epsg2proj(4326); %Proj4.0  projection string. , default is WGS 84 (corresponds to epsg 4326)
+
 	end
 	methods (Static)
@@ -30,5 +31,22 @@
 					self.shpfilename=getfieldvalue(options,'shpfilename','');
 					self.orientation=getfieldvalue(options,'orientation','normal');
-					self.epsg=getfieldvalue(options,'epsg',4326);
+
+					%figure out projection string: 
+					if exist(options,'epsg'),
+						if exist(options,'proj'),
+							error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
+						end
+						epsg=getfieldvalue(options,'epsg');
+						%convert into proj4: 
+						proj=epsg2proj(epsg);
+					elseif exist(options,'proj'),
+						if exist(options,'epsg'),
+							error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
+						end
+						proj=getfieldvalue(options,'proj');
+					else
+						proj= '+proj=longlat +datum=WGS84 +no_defs';
+					end
+					self.proj=proj;
 			end
 		end % }}}
@@ -37,5 +55,5 @@
 		self.shpfilename='';
 		self.orientation='normal';
-		self.epsg=4326;
+		self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
 		end % }}}
 		function disp(self) % {{{
@@ -45,5 +63,5 @@
 			fielddisplay(self,'shpfilename','shape file name');
 			fielddisplay(self,'orientation','orientation (default is ''normal'', can be ''reverse'')');
-			fielddisplay(self,'epsg','EPSG number defining projection for the shape file');
+			fielddisplay(self,'proj','shape file projection string (follows PROJ.4 standard)');
 
 		end % }}}
@@ -77,5 +95,4 @@
 			markersize=getfieldvalue(options,'markersize',1);
 			unitmultiplier=getfieldvalue(options,'unit',1);
-			epsg=getfieldvalue(options,'epsg',4326);
 			radius=getfieldvalue(options,'radius',6371012);
 			aboveground=getfieldvalue(options,'aboveground',1000);
@@ -84,4 +101,15 @@
 			label=getfieldvalue(options,'label','on');
 
+			if (exist(options,'epsg') &  exist(options,'proj')), 
+				error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
+			end
+			if exist(options,'epsg'),
+				proj=epsg2proj(getfieldvalue(options,'epsg'));
+			elseif exist(options,'proj'),
+				proj=getfieldvalue(options,'proj');
+			else
+				proj=epsg2proj(4326);
+			end
+
 			%read domain:
 			[path,name,ext]=fileparts(self.shpfilename);
@@ -92,8 +120,7 @@
 
 			%convert boundary to another reference frame:  {{{
-
 			for i=1:length(domain),
 				try, 
-					[x,y] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),sprintf('EPSG:%i',epsg));
+					[x,y] = gdaltransform(domain(i).x,domain(i).y,self.proj,proj);
 				catch me
 					disp(me.message());
@@ -125,4 +152,25 @@
 			%}}}
 		end % }}}
+		function checkconsistency(self,varargin) % {{{
+			%recover options
+			options=pairoptions(varargin{:});
+			tolerance=getfieldvalue(options,'tolerance',1e-5);
+		
+			%recover edges: 
+			edges=self.edges();
+
+			if strcmpi(edges.Geometry,'Point'),
+				return;
+			else
+				%check we don't have two identical vertices: 
+				x=edges.x; y=edges.y; 
+				distance=sqrt( (x(1:end-1)-x(2:end)).^2+(y(1:end-1)-y(2:end)).^2); 
+				dmax=max(distance);
+				isidentical=find(distance<dmax*tolerance);
+				if ~isempty(isidentical),
+					error(['boundary ' shpfilename ' has two  vertices extremely close to one another']);
+				end
+			end
+		end % }}}
 		function plot3d(self,varargin) % {{{
 			%recover options
@@ -136,5 +184,4 @@
 			markersize=getfieldvalue(options,'markersize',1);
 			unitmultiplier=getfieldvalue(options,'unit',1);
-			epsg=getfieldvalue(options,'epsg',4326);
 			radius=getfieldvalue(options,'radius',6371012);
 			aboveground=getfieldvalue(options,'aboveground',1000);
@@ -150,80 +197,38 @@
 			domain=shpread([self.shppath '/' name '.' ext]);
 
-			if epsg==4326,
-				%convert boundary to lat,long: {{{
-
-				for i=1:length(domain),
-					try, 
-						[lat,long] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),'EPSG:4326');
-					catch me
-						disp(me.message());
-						self.disp();
-					end
-					domain(i).x=long; domain(i).y=lat;
-				end
-
-				for i=1:length(domain),
-
-					%make sure lat,long are what they are supposed to be: 
-					%if any(domain(i).x>90 | domain(i).x<-90), 
-					%	long=domain(i).x; lat=domain(i).y;
-					%else
-					%	long=domain(i).y; lat=domain(i).x;
-					%end
-
-					%project on x,y,z reference frame.
-					[x,y,z]=AboveGround(domain(i).x,domain(i).y,radius,aboveground);
-					[xt,yt,zt]=AboveGround(domain(i).x+offset,domain(i).y+offset,radius,300000);
-					hold on;
-					if length(x)==1,
-						p=plot3(x,y,z,'k*'); 
-						set(p,'MarkerSize',markersize);
-						if strcmpi(label,'on'),
-							t=text(xt,yt,zt,self.shpfilename,'FontSize',fontsize);
-						end
-					else
-						p=plot3(x,y,z,'k-'); 
-						mid=floor(length(xt)/2);
-						if strcmpi(label,'on'),
-							text(xt(mid),yt(mid),zt(mid),self.shpfilename,'FontSize',fontsize);
-						end
-					end
-					set(p,'Color',color);
-					set(p,'LineWidth',linewidth);
-
-				end
-				%}}}
-			else
-				%convert boundary to another reference frame:  {{{
-
-				for i=1:length(domain),
-					try, 
-						[x,y] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),sprintf('EPSG:%i',epsg));
-					catch me
-						disp(me.message());
-						self.disp();
-					end
-					domain(i).x=x; domain(i).y=y;
-				end
-
-				for i=1:length(domain),
-					hold on;
-					if length(x)==1,
-						p=plot(x,y,'k*'); 
-						set(p,'MarkerSize',markersize);
-						if strcmpi(label,'on'),
-							t=text(x,y,self.shpfilename,'FontSize',fontsize);
-						end
-					else
-						p=plot(x,y,'k-'); 
-						if strcmpi(label,'on'),
-							text(sum(x)/length(x),sum(y)/length(y),self.shpfilename,'FontSize',fontsize);
-						end
-					end
-					set(p,'Color',color);
-					set(p,'LineWidth',linewidth);
-				end
-				%}}}
-			end
+			for i=1:length(domain),
+				try, 
+					[lat,long] = gdaltransform(domain(i).x,domain(i).y,self.proj,'EPSG:4326');
+				catch me
+					disp(me.message());
+					self.disp();
+				end
+				domain(i).x=long; domain(i).y=lat;
+			end
+
+			for i=1:length(domain),
+
+				%project on x,y,z reference frame.
+				[x,y,z]=AboveGround(domain(i).x,domain(i).y,radius,aboveground);
+				[xt,yt,zt]=AboveGround(domain(i).x+offset,domain(i).y+offset,radius,300000);
+				hold on;
+				if length(x)==1,
+					p=plot3(x,y,z,'k*'); 
+					set(p,'MarkerSize',markersize);
+					if strcmpi(label,'on'),
+						t=text(xt,yt,zt,self.shpfilename,'FontSize',fontsize);
+					end
+				else
+					p=plot3(x,y,z,'k-'); 
+					mid=floor(length(xt)/2);
+					if strcmpi(label,'on'),
+						text(xt(mid),yt(mid),zt(mid),self.shpfilename,'FontSize',fontsize);
+					end
+				end
+				set(p,'Color',color);
+				set(p,'LineWidth',linewidth);
+
+			end
+
 		end % }}}
 	end
Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 24778)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 24779)
@@ -23,4 +23,5 @@
 		transitions      = {};
 		eltransitions      = {};
+		planet           = '';
 		%}}}
 	end
@@ -39,5 +40,9 @@
 				
 				%recover the earth model:
-				slm.earth = getfieldvalue(options,'earth');
+				slm.earth = getfieldvalue(options,'earth',0);
+
+				%set planet type: 
+				slm.planet=getfieldvalue(options,'planet','earth');
+
 			end
 		end
@@ -91,5 +96,6 @@
 			slm.range             = {};
 			slm.transitions       = {};
-			slm.eltransitions       = {};
+			slm.eltransitions     = {};
+			slm.planet            = 'earth';
 		end
 		%}}}
@@ -153,4 +159,7 @@
 		options=pairoptions(varargin{:}); 
 		force=getfieldvalue(options,'force',0);
+		
+		%initialize, to avoid issues of having more transitions than meshes. 
+		self.transitions={}; self.eltransitions={};
 
 		%for elements:
@@ -161,5 +170,5 @@
 		for i=1:length(self.icecaps),
 			mdi=self.icecaps{i};
-			mdi=TwoDToThreeD(mdi);
+			mdi=TwoDToThreeD(mdi,self.planet);
 	
 			%for elements: 
@@ -182,4 +191,10 @@
 		end
 		plotmodel(self.earth,'data',flags,'coastline','on');
+
+		end % }}}
+		function checkbasinconsistency(self) % {{{
+			for i=1:length(self.basins),
+				self.basins{i}.checkconsistency();
+			end
 
 		end % }}}
@@ -300,5 +315,5 @@
 			models=self.icecaps;
 			for i=1:length(models),
-				models{i}=TwoDToThreeD(models{i});
+				models{i}=TwoDToThreeD(models{i},self.planet);
 			end
 			
@@ -350,5 +365,5 @@
 		function transfer(self,string) % {{{
 			%Recover field size in one icecap: 
-			eval(['n=length(self.icecaps{1}.' string ');']);
+			eval(['n=size(self.icecaps{1}.' string ',1);']);
 			if n==self.icecaps{1}.mesh.numberofvertices,
 				eval(['self.earth.' string '=zeros(self.earth.mesh.numberofvertices,1);']);
@@ -356,4 +371,33 @@
 					eval(['self.earth.' string '(self.transitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
 				end
+			elseif n==(self.icecaps{1}.mesh.numberofvertices+1),
+				%dealing with a transient dataset. 
+				%check that all timetags are similar between all icecaps:  %{{{
+				for i=1:length(self.icecaps),
+					eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
+					for j=(i+1):length(self.icecaps),
+						eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
+						if ~isequal(capfieldi(end,:),capfieldj(end,:)),
+							error(['Time stamps for ' string ' field is different between icecaps ' num2str(i) ' and ' num2str(j)]);
+						end
+					end
+				end
+				eval(['capfield1= self.icecaps{1}.' string ';']); 
+				times=capfield1(end,:);
+				nsteps=length(times);
+				%}}}
+				%initialize:  %{{{
+				eval(['field=zeros(self.earth.mesh.numberofvertices+1,' num2str(nsteps) ');']);
+				eval(['field(end,:)=times;']); %transfer the times only, not the values
+				%}}}
+				%transfer all time fields: {{{
+				for i=1:length(self.icecaps),
+					eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
+					for j=1:nsteps,
+						eval(['field(self.transitions{' num2str(i) '},' num2str(j) ')=capfieldi(1:end-1,' num2str(j) ');']); %transfer only the values, not the time.
+					end
+				end
+				eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location
+				%}}}
 			elseif n==self.icecaps{1} .mesh.numberofelements,
 				eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']);
Index: /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m	(revision 24778)
+++ /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m	(revision 24779)
@@ -1,11 +1,11 @@
-function md=TwoDToThreeD(md)
+function md=TwoDToThreeD(md,planet)
 
 	%reproject model into lat,long if necessary:
-	if md.mesh.epsg~=4326,
-		[md.mesh.x,md.mesh.y]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),sprintf('EPSG:4326'));
+	if ~strcmpi(md.mesh.proj,epsg2proj(4326)),
+		[md.mesh.x,md.mesh.y]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
 	end
 
 	%Make a 3dsurface mesh out of this: 
-	R=6.371012*10^6; 
+	R=planetradius(planet);
 
 	%we assume x and y hold the long,lat values:
@@ -14,4 +14,5 @@
 	mask=md.mask;
 
+	%Assume spherical body: 
 	x = R .* cosd(lat) .* cosd(long);
 	y = R .* cosd(lat) .* sind(long);
Index: /issm/trunk-jpl/src/m/mesh/meshintersect3d.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshintersect3d.m	(revision 24778)
+++ /issm/trunk-jpl/src/m/mesh/meshintersect3d.m	(revision 24779)
@@ -66,4 +66,6 @@
 	end
 
+	if(~isempty(find(indices==0))) error('issue with transition vector having one empty slot'); end;
+
 
 %		if length(s)>1,