Index: /issm/trunk/src/m/classes/public/BasinConstrain.m
===================================================================
--- /issm/trunk/src/m/classes/public/BasinConstrain.m	(revision 1337)
+++ /issm/trunk/src/m/classes/public/BasinConstrain.m	(revision 1338)
@@ -61,12 +61,4 @@
 md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
 
-
-%make sure any grid with NaN velocity is spc'd:
-pos=find(isnan(md.vel_obs_raw));
-md.gridondirichlet_diag(pos)=1;
-%we spc to the smoothed value, so that control methods don't go berserk trying to figure out what reaction force to apply for the spc to stand.
-md.dirichletvalues_diag(pos,1)=md.vx_obs(pos); 
-md.dirichletvalues_diag(pos,2)=md.vy_obs(pos); 
-
 %make sure icefronts that are completely spc'd are taken out:
 free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);
Index: /issm/trunk/src/m/classes/public/extrude.m
===================================================================
--- /issm/trunk/src/m/classes/public/extrude.m	(revision 1337)
+++ /issm/trunk/src/m/classes/public/extrude.m	(revision 1338)
@@ -36,7 +36,7 @@
 end
 
-if extrusionexponent<0,
+if extrusionexponent<=0,
 	help extrude;
-	error('extrusionexponent must be >0');
+	error('extrusionexponent must be >=0');
 end
 
Index: /issm/trunk/src/m/classes/public/mesh/meshyams.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1337)
+++ /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1338)
@@ -1,49 +1,48 @@
-function md=meshyams(md,velpath,domainoutline,varargin);
+function md=meshyams(md,varargin);
 %MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
 %
 %   Usage:
-%      md=meshyams;
+%      md=meshyams(md,varargin);
+%      where varargin is a lit of paired arguments. 
+%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
+%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
+%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
+%                          this option is used to minimize the metric on water (no refinement)
+%      optional arguments: 'resolution': initial mesh resolution [m]
+%      optional arguments: 'nsteps': number of steps of mesh adaptation
+%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
+%      optional arguments: 'hmin': minimum edge length
+%      optional arguments: 'hmanx': maximum edge
+%      
+%
+%   Examples:
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat');
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
 
-%PIG settings
-%{
-nsteps=4;
-resolution=10000;
-hmin=1500;
-hmax=10^7;
-gradation=3*ones(nsteps,1);
-epsilon=10^-0;
-scale=1;
-%}
+%recover options
+yamsoptions=recover_options(varargin{:});
 
-%%{
-%Antarctica settings
-nsteps=6;
-resolution=5000;
-hmin=300;        %300m
-hmax=150*10^3;   %150km
-gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
-epsilon=2*10^-0; %3m/a interpolation error
-scale=1;
-%%}
+%add default options
+yamsoptions=process_mesh_options(yamsoptions);
 
-%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
+%recover some fields
+domainoutline=yamsoptions.domainoutline;
+groundeddomain=yamsoptions.groundeddomain;
+velocities=yamsoptions.velocities;
+resolution=yamsoptions.resolution;
+gradation=yamsoptions.gradation;
+nsteps=yamsoptions.nsteps;
+epsilon=yamsoptions.epsilon;
+hmin=yamsoptions.hmin;
+hmax=yamsoptions.hmax;
 
 %mesh with initial resolution
+disp('Initial mesh generation...');
 md=mesh(md,domainoutline,resolution);
 
 %load velocities 
-load(velpath);
-
-disp(['First mesh, number of elements: ' num2str(md.numberofelements)]);
+disp('loading velocities...');
+load(velocities);
 
 %start mesh adaptation
@@ -53,11 +52,11 @@
 	%interpolate velocities onto mesh
 	disp('   interpolating velocities');
-	md.vx_obs=InterpFromGrid(x_m,y_m,vx,md.x,md.y,0)*scale;
-	md.vy_obs=InterpFromGrid(x_m,y_m,vy,md.x,md.y,0)*scale;
-	md.vel_obs=averaging(md,sqrt(md.vx_obs.^2+md.vy_obs.^2),2);
+	vx_obs=InterpFromGrid(x_m,y_m,vx,md.x,md.y,0);
+	vy_obs=InterpFromGrid(x_m,y_m,vy,md.x,md.y,0);
+	field=averaging(md,sqrt(vx_obs.^2+vy_obs.^2),2);
 
 	%set gridonwater field
-	if waterflag,
-		gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundedoutline,1),'node',2);
+	if ~strcmp(groundeddomain,'N/A'),
+		gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundeddomain,1),'node',2);
 		md.gridonwater=ones(md.numberofgrids,1);
 		md.gridonwater(find(gridground))=0;
@@ -68,5 +67,5 @@
 	%adapt according to velocities
 	disp('   adapting');
-	md=YamsCall(md,md.vel_obs,hmin,hmax,gradation(i),epsilon);
+	md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
 end
 	
@@ -78,4 +77,5 @@
 
 %Recreate the segments
+disp('Creating segments...');
 elementconnectivity=md.elementconnectivity;
 elementonboundary=double(elementconnectivity(:,end)==0);
Index: /issm/trunk/src/m/classes/public/mesh/process_mesh_options.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/process_mesh_options.m	(revision 1338)
+++ /issm/trunk/src/m/classes/public/mesh/process_mesh_options.m	(revision 1338)
@@ -0,0 +1,125 @@
+function outoptions=process_mesh_options(options)
+%DEFAULT_MESH_OPTIONS - set up default options for mesh phase
+%
+%   Usage:
+%      options=process_mesh_options(options)
+%
+%   See also: MESH,MESHYAMS,RECOVER_OPTIONS
+
+%domainoutline: check on this option, error out otherwise
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'domainoutline'),
+		if exist(options{i,2}),
+			domainoutline=options{i,2};
+			found=1;
+		else
+			error(['process_mesh_options error message: domain outline file ' options{i,2} ' not found!']);
+		end
+	end
+end
+if ~found,
+	error('process_mesh_options error message: no ''domainoutline'' was provided');
+end
+
+%velocities: check on this option, error out otherwise
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'velocities'),
+		if exist(options{i,2}),
+			velocities=options{i,2};
+			found=1;
+		else
+			error(['process_mesh_options error message: velocities file ' options{i,2} ' not found!']);
+		end
+	end
+end
+if ~found,
+	error('process_mesh_options error message: no ''velocities'' was provided');
+end
+
+%groundeddomain: check on this option, error out otherwise
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'groundeddomain'),
+		if exist(options{i,2}),
+			groundeddomain=options{i,2};
+			found=1;
+		else
+			error(['process_mesh_options error message: grounded domain file ' options{i,2} ' not found!']);
+		end
+	end
+end
+if ~found,
+	groundeddomain='N/A';
+end
+
+%nsteps
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'nsteps'),
+		nsteps=options{i,2};
+		found=1;
+	end
+end
+if ~found, nsteps=6;end
+
+%resolution
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'resolution'),
+		resolution=options{i,2};
+		found=1;
+	end
+end
+if ~found, resolution=5000;end
+
+%hmin
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'hmin'),
+		hmin=options{i,2};
+		found=1;
+	end
+end
+if ~found, hmin=300;end
+
+%hmax
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'hmax'),
+		hmax=options{i,2};
+		found=1;
+	end
+end
+if ~found, hmax=150*10^3;end
+
+%epsilon
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'epsilon'),
+		epsilon=options{i,2};
+		found=1;
+	end
+end
+if ~found, epsilon=2;end
+
+%gradation
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'gradation'),
+		gradation=options{i,2};
+		found=1;
+	end
+end
+if ~found, 
+	gradation=[1.5*ones(2,1);3*ones(nsteps,1)];
+	gradation(nsteps+1:end)=[];
+end
+
+if length(gradation)~=nsteps
+	error(['process_mesh_options error message: gradation should have the same length as nsteps (' num2str(nsteps) ')']);
+end
+
+%output options
+outoptions=struct('domainoutline',domainoutline,'velocities',velocities,'groundeddomain',groundeddomain,'nsteps',nsteps,'resolution',resolution,'hmin',hmin,'hmax',hmax,'epsilon',epsilon,'gradation',gradation);
Index: /issm/trunk/src/m/classes/public/mesh/recover_mesh_options.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/recover_mesh_options.m	(revision 1338)
+++ /issm/trunk/src/m/classes/public/mesh/recover_mesh_options.m	(revision 1338)
@@ -0,0 +1,22 @@
+function options=recover_mesh_options(md,varargin)
+%RECOVER_mesh_OPTIONS - recover solution options
+%
+%   Usage:
+%      options=recover_mesh_options(md,varargin);
+%
+%   See also: mesh
+
+%initialize options.
+options=cell(0,2);
+
+%make sure length(varargin) is even, ie options come in pairs.
+if mod(length(varargin),2),
+	error('recover_mesh_options error message: an even number of options is necessary');
+end
+
+%go through varargin, extract options 
+for i=1:length(varargin)/2,
+
+	options(end+1,:)={varargin{2*i-1} varargin{2*i}};
+
+end
Index: /issm/trunk/src/m/classes/public/recover_areas.m
===================================================================
--- /issm/trunk/src/m/classes/public/recover_areas.m	(revision 1337)
+++ /issm/trunk/src/m/classes/public/recover_areas.m	(revision 1338)
@@ -15,5 +15,4 @@
 
 	%make sure length(varargin) is even
-	varargin=varargin{1};
 	if mod(length(varargin),2),
 		error('recover_elementstype_options error message: an even number of options is necessary');
Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 1337)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 1338)
@@ -26,5 +26,5 @@
 end
 
-[hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
+[hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin{:});
 
 %Flag the elements that has not been flagged as filltype
Index: /issm/trunk/src/m/classes/public/solve.m
===================================================================
--- /issm/trunk/src/m/classes/public/solve.m	(revision 1337)
+++ /issm/trunk/src/m/classes/public/solve.m	(revision 1338)
@@ -20,5 +20,5 @@
 
 %recover options
-options=recover_solve_options(md,varargin{:});
+options=recover_options(varargin{:});
 
 %add default options
Index: /issm/trunk/src/m/utils/Miscellaneous/recover_options.m
===================================================================
--- /issm/trunk/src/m/utils/Miscellaneous/recover_options.m	(revision 1338)
+++ /issm/trunk/src/m/utils/Miscellaneous/recover_options.m	(revision 1338)
@@ -0,0 +1,22 @@
+function options=recover_options(varargin)
+%RECOVER_OPTIONS - recover options from a list of inputs
+%
+%   Usage:
+%      options=recover_solve_options(varargin);
+%
+%   See also: RECOVER_PLOT_OPTIONS, MESHYAMS,SOLVE
+
+%initialize options.
+options=cell(0,2);
+
+%make sure length(varargin) is even, ie options come in pairs.
+if mod(length(varargin),2),
+	error('recover_options error message: an even number of options is necessary');
+end
+
+%go through varargin, extract options 
+for i=1:length(varargin)/2,
+
+	options(end+1,:)={varargin{2*i-1} varargin{2*i}};
+
+end