Index: /issm/trunk/src/m/classes/public/mesh/meshyams.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 3259)
+++ /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 3260)
@@ -55,5 +55,5 @@
 else
 	md=mesh(md,domainoutline,riftoutline,resolution);
-	md=meshprocessrifts(md);
+	md=meshprocessrifts(md,domainoutline);
 end
 disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);
Index: /issm/trunk/src/m/classes/public/mesh/rifts/meshprocessrifts.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/meshprocessrifts.m	(revision 3259)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/meshprocessrifts.m	(revision 3260)
@@ -1,16 +1,26 @@
-function md=meshprocessrifts(md)
+function md=meshprocessrifts(md,domainoutline)
 %MESHPROCESSRIFTS - process mesh when rifts are present
 %
 %   split rifts inside mesh (rifts are defined by presence of
 %   segments inside the domain outline)
+%   if domain outline is provided, check for rifts that could touch it, and open them up.
 %
 %   Usage:
 %      md=meshprocessrifts(md)
+%
+%   Ex: 
+%      md=meshprocessrifts(md);
+%      md=meshprocessrifts(md,'DomainOutline.exp');
+%
 
-%some checks on list of arguments
-if ((nargin~=1) | (nargout~=1)),
-	help meshprocessrifts
-	error('meshprocessrifts error messagei: bad usage');
+%some checks on arguments: 
+if nargout~=1,
+	error('meshprocessrifts usage error:');
 end
+
+if nargin~=2,
+	error('meshprocessrifts usage error:');
+end
+
 
 [md.elements,md.x,md.y,md.segments,md.segmentmarkers,md.rifts]=TriMeshProcessRifts(md.elements,md.x,md.y,md.segments,md.segmentmarkers);
@@ -33,2 +43,24 @@
 	md.rifts(i).tip2coordinates=[md.x(md.rifts(i).tips(2)) md.y(md.rifts(i).tips(2))];
 end
+
+%In case we have rifts that open up the domain outline, we need to open them: 
+flags=ContourToMesh(md.elements,md.x,md.y,expread(domainoutline,1),'node',0);
+found=0;
+for i=1:md.numrifts,
+	if flags(md.rifts(i).tips(1))==0,
+		found=1;
+		break;
+	end
+	if flags(md.rifts(i).tips(2))==0,
+		found=1;
+		break;
+	end
+end
+if found,
+	md=meshprocessoutsiderifts(md,domainoutline);
+end
+
+%get elements that are not correctly oriented in the correct direction:
+aires=GetAreas(md.elements,md.x,md.y);
+pos=find(aires<0);
+md.elements(pos,:)=[md.elements(pos,2) md.elements(pos,1) md.elements(pos,3)];
Index: /issm/trunk/src/m/classes/public/mesh/rifts/meshyamsrecreateriftsegments.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/meshyamsrecreateriftsegments.m	(revision 3259)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/meshyamsrecreateriftsegments.m	(revision 3260)
@@ -7,26 +7,82 @@
 			rift=md.rifts(i);
 
-			%find tip1 and tip2 for this rift, in the new mesh created by yams.
-			pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
-			tip1=md.segments(pos,1);
-			pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
-			tip2=md.segments(pos,1);
+			%closed rifts first:
+			if length(rift.tips)==2,
 
-			%start from tip1, and build segments of this rift. 
-			pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
-			pos_record=[pos_record; pos];
-			riftsegs=md.segments(pos,:);
-			while 1,
-				A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
-				%find other segment that holds B.
-				pos=find(md.segments(:,1)==B);
+				%find tip1 and tip2 for this rift, in the new mesh created by yams.
+				pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+				tip1=md.segments(pos,1);
+				pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
+				tip2=md.segments(pos,1);
+
+				%start from tip1, and build segments of this rift. 
+				pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+				pos_record=[pos_record; pos];
+				riftsegs=md.segments(pos,:);
+				while 1,
+					A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
+					%find other segment that holds B.
+					pos=find(md.segments(:,1)==B);
+					pos_record=[pos_record; pos];
+					riftsegs=[riftsegs; md.segments(pos,:)];
+					if riftsegs(end,2)==tip1, 
+						break;
+					end
+				end
+				md.rifts(i).segments=riftsegs;
+				md.rifts(i).tips=[tip1 tip2];
+
+			else
+				%ok, this is a rift that opens up to the domain outline.  One tip is going to be 
+				%double, the other one, single. We are going to start from the single tip, towards the two 
+				%other doubles
+
+				%find tip1 and tip2 for this rift, in the new mesh created by yams.
+				pos1=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+				tip1=md.segments(pos1,1);
+				pos2=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
+				tip2=md.segments(pos2,1);
+				if length(tip1)==2,
+					%swap.
+					temp=tip1; tip1=tip2; tip2=temp;
+					temp=pos1; pos1=pos2; pos2=temp;
+					pos=pos1;
+				else
+					pos=pos1;
+				end
+
+				pos_record=[pos_record; pos];
+				riftsegs=md.segments(pos,:);
+				while 1,
+					A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
+					%find other segment that holds B.
+					pos=find(md.segments(:,1)==B);
+					pos_record=[pos_record; pos];
+					riftsegs=[riftsegs; md.segments(pos,:)];
+					if ((riftsegs(end,2)==tip2(1)) | (riftsegs(end,2)==tip2(2))), 
+						%figure out which tip we reached
+						if riftsegs(end,2)==tip2(1), index=2; else index=1; end
+						break;
+					end
+				end
+
+				%ok, now, we start from the other tip2, towards tip1
+				pos=pos2(index);
 				pos_record=[pos_record; pos];
 				riftsegs=[riftsegs; md.segments(pos,:)];
-				if riftsegs(end,2)==tip1, 
-					break;
+				while 1,
+					A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
+					%find other segment that holds B.
+					pos=find(md.segments(:,1)==B);
+					pos_record=[pos_record; pos];
+					riftsegs=[riftsegs; md.segments(pos,:)];
+					if riftsegs(end,2)==tip1, 
+						break;
+					end
 				end
+				md.rifts(i).segments=riftsegs;
+				md.rifts(i).tips=[tip1 tip2(1) tip2(2)];
+
 			end
-			md.rifts(i).segments=riftsegs;
-			md.rifts(i).tips=[tip1 tip2];
 		end
 	end
Index: /issm/trunk/src/m/classes/public/modis.m
===================================================================
--- /issm/trunk/src/m/classes/public/modis.m	(revision 3260)
+++ /issm/trunk/src/m/classes/public/modis.m	(revision 3260)
@@ -0,0 +1,37 @@
+function [xm,ym,modis]=modis(modisgeotif,xlim,ylim)
+%MODIS - from modis geotiff, return image
+%
+%   Usage:
+%      [xm,ym,modis]=modis(modisgeotif,xlim,ylim)
+%
+
+global ISSM_DIR
+
+%find gdal coordinates
+x0=min(xlim);
+x1=max(xlim);
+
+y0=min(ylim);
+y1=max(ylim);
+
+%Get path  to gdal binaries
+path_gdal=[ISSM_DIR '/externalpackages/gdal/install/bin/'];
+
+%Was gdal compiled? 
+if ~exist([path_gdal 'gdal_translate']),
+	error(['modis error message: GDAL library needs to be compiled to use this routine. Compile GDAL in ' ISSM_DIR '/externalpackages/gdal to use this routine.']);
+end
+
+inputname='./temp.tif';
+system([path_gdal 'gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' modisgeotif ' ' inputname ]);
+
+%Read in temp.tif:
+modis=double(flipud(imread('temp.tif','TIFF')));
+xm=(x0:(x1-x0)/(size(md.sarpwr,2)-1):x1);
+ym=(y0:(y1-y0)/(size(md.sarpwr,1)-1):y1);
+
+%Erase image
+system('rm -rf ./temp.tif');
+
+
+end
Index: /issm/trunk/src/m/classes/public/plot/plot_rifts.m
===================================================================
--- /issm/trunk/src/m/classes/public/plot/plot_rifts.m	(revision 3259)
+++ /issm/trunk/src/m/classes/public/plot/plot_rifts.m	(revision 3260)
@@ -17,5 +17,5 @@
 if isstruct(md.rifts),
 	
-	for i=1:size(md.rifts,1),
+	for i=1:md.numrifts,
 		penaltypairs=md.rifts(i).penaltypairs;
 
@@ -26,4 +26,33 @@
 			x(penaltypairs(j,1))=x(penaltypairs(j,1))-normal(1)*offset;
 			y(penaltypairs(j,1))=y(penaltypairs(j,1))-normal(2)*offset;
+		end
+		if length(md.rifts(i).tips)==3,
+			tip=md.rifts(i).tips(3);
+			%who is tip connected to? 
+			if isconnected(md.elements,penaltypairs(1,1),tip),
+				normal(1)=penaltypairs(1,5);
+				normal(2)=penaltypairs(1,6);
+				x(tip)=x(tip)-normal(1)*offset;
+				y(tip)=y(tip)-normal(2)*offset;
+			end
+
+			if isconnected(md.elements,penaltypairs(1,2),tip),
+				normal(1)=penaltypairs(1,5);
+				normal(2)=penaltypairs(1,6);
+				x(tip)=x(tip)+normal(1)*offset;
+				y(tip)=y(tip)+normal(2)*offset;
+			end
+			if isconnected(md.elements,penaltypairs(end,1),tip),
+				normal(1)=penaltypairs(end,5);
+				normal(2)=penaltypairs(end,6);
+				x(tip)=x(tip)-normal(1)*offset;
+				y(tip)=y(tip)-normal(2)*offset;
+			end
+			if isconnected(md.elements,penaltypairs(end,2),tip),
+				normal(1)=penaltypairs(end,5);
+				normal(2)=penaltypairs(end,6);
+				x(tip)=x(tip)+normal(1)*offset;
+				y(tip)=y(tip)+normal(2)*offset;
+			end
 		end
 	end