Index: /issm/trunk/src/m/classes/public/mesh/meshbamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshbamg.m	(revision 3198)
+++ /issm/trunk/src/m/classes/public/mesh/meshbamg.m	(revision 3199)
@@ -6,10 +6,11 @@
 %      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]
+%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
+%      arguments can be: 'thickness': matlab file containing the thickness [m]
 %      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: 'epsilon': average interpolation error wished [m/yr [m]]
 %      optional arguments: 'hmin': minimum edge length
 %      optional arguments: 'hmanx': maximum edge
@@ -20,4 +21,5 @@
 %      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);
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','thickness','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
 
 %recover options
@@ -26,5 +28,5 @@
 
 %recover some fields
-disp('MeshYams Options:')
+disp('MeshBamg Options:')
 domainoutline=getfieldvalueerr(options,'domainoutline');
 disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
@@ -33,4 +35,6 @@
 velocities=getfieldvalueerr(options,'velocities');
 disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
+thickness=getfieldvalue(options,'thickness','none');
+disp(sprintf('   %-15s: ''%s''','Thickness',thickness));
 resolution=getfieldvalue(options,'resolution',5000);
 disp(sprintf('   %-15s: %f','Resolution',resolution));
@@ -53,6 +57,16 @@
 %load velocities 
 disp('loading velocities...');
-Names=VelFindVarNames(velocities);
+NamesV=VelFindVarNames(velocities);
 Vel=load(velocities);
+
+%compute thickness?
+if ~strcmpi(thickness,'none'),
+	disp('loading thickness...');
+	NamesH=FieldFindVarNames(thickness);
+	Thi=load(thickness);
+	thicknesspresent=1;
+else
+	thicknesspresent=0;
+end
 
 %start mesh adaptation
@@ -62,13 +76,23 @@
 	%interpolate velocities onto mesh
 	disp('   interpolating velocities...');
-	if strcmpi(Names.interp,'grid'),
-		vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
-		vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
+	if strcmpi(NamesV.interp,'grid'),
+		vx_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.x,md.y,0);
+		vy_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md.x,md.y,0);
 	else
-		vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
-		vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
+		vx_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.x,md.y,0);
+		vy_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md.x,md.y,0);
 	end
 	field=sqrt(vx_obs.^2+vy_obs.^2);
-
+	
+	if thicknesspresent,
+		disp('   interpolating thickness...');
+		if strcmpi(NamesH.interp,'grid'),
+			h=InterpFromGridToMesh(Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.x,md.y,0);
+		else
+			h=InterpFromMeshToMesh2d(Thi.(NamesH.indexname),Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.x,md.y,0);
+		end
+		field=[field h];
+	end
+	
 	%set gridonwater field
 	if ~strcmp(groundeddomain,'N/A'),