Index: /issm/trunk/src/m/utils/Mesh/MeshYams.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/MeshYams.m	(revision 1335)
+++ /issm/trunk/src/m/utils/Mesh/MeshYams.m	(revision 1336)
@@ -1,3 +1,3 @@
-function md=MeshYams(md,velpath,domainoutline);
+function md=MeshYams(md,velpath,domainoutline,varargin);
 %MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
 %
@@ -23,7 +23,19 @@
 hmax=150*10^3;   %150km
 gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
-epsilon=3*10^-0; %3m/a interpolation error
+epsilon=2*10^-0; %3m/a interpolation error
 scale=1;
 %%}
+
+%check number of arguments
+waterflag=0;
+if nargin==4,
+	groundedoutline=varargin{1};
+	if exist(groundedoutline);
+		disp(['grounded ice domain found. Metric will be minimum on water']);
+		waterflag=1;
+	else
+		error(['MeshYams error message: file ' groundedoutline ' not found.']);
+	end
+end
 
 %mesh with initial resolution
@@ -45,4 +57,13 @@
 	md.vel_obs=averaging(md,sqrt(md.vx_obs.^2+md.vy_obs.^2),2);
 
+	%set gridonwater field
+	if waterflag,
+		gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundedoutline,1),'node',2);
+		md.gridonwater=ones(md.numberofgrids,1);
+		md.gridonwater(find(gridground))=0;
+	else
+		md.gridonwater=zeros(md.numberofgrids,1);
+	end
+
 	%adapt according to velocities
 	disp('   adapting');
@@ -51,2 +72,41 @@
 	
 disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);
+
+%Now, build the connectivity tables for this mesh.
+md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
+
+%Recreate the segments
+elementconnectivity=md.elementconnectivity;
+elementonboundary=double(elementconnectivity(:,end)==0);
+pos=find(elementonboundary);
+num_segments=length(pos);
+segments=zeros(num_segments,3);
+for i=1:num_segments,
+	el1=pos(i);
+	els2=elementconnectivity(el1,find(elementconnectivity(el1,:)));
+	flag=intersect(md.elements(els2(1),:),md.elements(els2(2),:));
+	nods1=md.elements(el1,:);
+	nods1(find(nods1==flag))=[];
+	segments(i,:)=[nods1 el1];
+
+	ord1=find(nods1(1)==md.elements(el1,:));
+	ord2=find(nods1(2)==md.elements(el1,:));
+
+	%swap segment grids if necessary
+	if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
+		temp=segments(i,1);
+		segments(i,1)=segments(i,2);
+		segments(i,2)=temp;
+	end
+	segments(i,1:2)=fliplr(segments(i,1:2));
+end
+md.segments=segments;
+
+%Fill in rest of fields:
+md.z=zeros(md.numberofgrids,1);
+md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
+md.gridonbed=ones(md.numberofgrids,1);
+md.gridonsurface=ones(md.numberofgrids,1);
+md.elementonbed=ones(md.numberofelements,1);
+md.elementonsurface=ones(md.numberofelements,1);
Index: /issm/trunk/src/m/utils/Mesh/YamsCall.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/YamsCall.m	(revision 1335)
+++ /issm/trunk/src/m/utils/Mesh/YamsCall.m	(revision 1336)
@@ -88,6 +88,12 @@
 
 %some corrections for 0 eigen values
-metric(pos1,:)=repmat([1/hmin^2 0 1/hmin^2],length(pos1),1);
-metric(pos2,:)=repmat([1/hmin^2 0 1/hmin^2],length(pos2),1);
+metric(pos1,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos1),1);
+metric(pos2,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos2),1);
+
+%take care of water elements
+if length(md.gridonwater)==numberofgrids;
+	pos=find(md.gridonwater);
+	metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
+end
 
 if any(isnan(metric)),
@@ -143,5 +149,5 @@
 md.numberofgrids=size(Coor,1);
 md.numberofelements=size(Tria,1);
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
 %clean up: