function md=BuildAdaptedMesh
%BUILDADAPTEDMESH - Build model of Antarctica by refining according to observed velocity error estimator
%
%   Usage:
%      md=BuildAdaptedMesh

nsteps=1;
resolution=50000;
scaling=1.2;
epsilon=.5;

velpath='/proj/ice/larour/Glaciology/Model_Data/RignotAntarcticaVel1km/RignotAntVel';
%velpath='/proj/ice/larour/Glaciology/Model_Data/BamberAntarcticaVel1km/velocity';
%velpath='/proj/ice/larour/Glaciology/Model_Data/mosaicAnt1km_ramp_ers_alos/RignotAntVel';

domainoutline='DomainOutline.exp';
groundingline='GroundingLine.exp';
icestreams='Icestream.exp';
parameterfile='Basins.par';

%clear data
clear md

%build new model
md=model;

%mesh with initial resolution
md=mesh(md,domainoutline,resolution);
disp(['First mesh, number of elements: ' num2str(md.numberofelements)]);

%load velocities 
eval(['load ' velpath]);

%start mesh adaptation
for i=1:nsteps,
	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);

	%interpolate velocities onto mesh
	disp('   interpolating velocities');
	md.vx_obs=DataInterp(x_m,y_m,vx*md.yts,md.x,md.y);
	md.vy_obs=DataInterp(x_m,y_m,vy*md.yts,md.x,md.y);
	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);

	%adapt according to velocities
	disp('   adapting');
	md=meshadaptation(md,md.vel_obs,scaling,epsilon);
end
	
disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);

md.vx_obs=DataInterp(x_m,y_m,vx*md.yts,md.x,md.y);
md.vy_obs=DataInterp(x_m,y_m,vy*md.yts,md.x,md.y);
md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
return;

%parameterize model
disp('Parameterizing');
md=geography(md,groundingline,icestreams);
md=parameterize(md,parameterfile);

%md.vx_obs=DataInterp(x_m,y_m,vx*md.yts,md.x,md.y);
%md.vy_obs=DataInterp(x_m,y_m,vy*md.yts,md.x,md.y);
%md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
%return;
