Index: /issm/trunk/src/m/classes/public/geography2.m
===================================================================
--- /issm/trunk/src/m/classes/public/geography2.m	(revision 2985)
+++ /issm/trunk/src/m/classes/public/geography2.m	(revision 2986)
@@ -1,3 +1,10 @@
 function md=geography2(md,landname,iceshelfname,icesheetname)
+%GEOGRAPHY2 - establish land, ice sheet and ice shelf areas in a domains.
+%
+%   Usage:
+%      md=geography2(md,landname,iceshelfname,icesheetname)
+%
+%   Examples:
+%      md=geography2(md,'LandName.exp','Iceshelves.exp','Islands.exp');
 
 %Get assigned fields
@@ -7,5 +14,16 @@
 
 %recover elements and grids on land.
-[gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
+if ischar(landname),
+	[gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
+elseif isfloat(landname),
+	if size(landname,1)~=md.numberofelements,
+		error('Landname for area must be of same size as number of elements in model');
+	end
+	elementonland=landname;
+	gridonland=zeros(md.numberofgrids,1);
+	gridonland(md.elements(find(elementonland),:))=1;
+else
+	error('Invalid area option option');
+end
 
 %any element with 3 grids on land should be on land:
@@ -42,21 +60,15 @@
 
 %recover arrays of ice shelf grids and elements, and ice sheet grids and elements.
-if strcmp(iceshelfname,''), %no iceshelf contour file, we are dealing with a pure ice sheet.
-	gridoniceshelf=zeros(md.numberofgrids,1);
-	elementoniceshelf=zeros(md.numberofelements,1);
-elseif strcmp(iceshelfname,'all'), %we are dealing with a pure ice shelf.
-	gridoniceshelf=ones(md.numberofgrids,1);
-	elementoniceshelf=ones(md.numberofelements,1);
-else
-	[gridoniceshelf,elementoniceshelf]=ContourToMesh(elements,x,y,expread(iceshelfname,1),'element and node',2);
-end
+elementoniceshelf=FlagElements(md,iceshelfname);
+elementonicesheet=FlagElements(md,icesheetname);
 
-if strcmp(icesheetname,''), %no icesheet contour file, we are dealing with a pure ice shelf.
-	gridonicesheet=zeros(md.numberofgrids,1);
-	elementonicesheet=zeros(md.numberofelements,1);
-else
-	[gridonicesheet,elementonicesheet]=ContourToMesh(elements,x,y,expread(icesheetname,1),'element and node',2);
-end
-
+%Because icesheet grids and elements can be included into an iceshelf, we need to update. Remember, all the previous 
+%arrays come from domain outlines that can intersect one another: 
+gridoniceshelf=zeros(md.numberofgrids,1);
+gridonicesheet=zeros(md.numberofgrids,1);
+elementoniceshelf=double((elementoniceshelf & ~elementonicesheet));
+elementonicesheet=double(~elementoniceshelf);
+gridoniceshelf(md.elements(find(elementoniceshelf),:))=1;
+gridonicesheet(md.elements(find(elementonicesheet),:))=1;
 
 %now correct, so that none of the iceshelf and icesheet elements and nodes are in the water.