steps=[1 2 3 4];

%Hard coded parameters
cluster=generic('name',oshostname,'np',2);
ncdata='../../Data/Greenland_5km_dev1.2.nc';

%Run Steps
if any(steps==1)

	disp('	Step 1: Mesh creation'); 
	md=triangle(model,'Jks3_trimmed.exp',2000);

	%Get observed velocity field on mesh nodes
	x1=double(ncread(ncdata,'x1'));
	y1=double(ncread(ncdata,'y1'));
	velx=double(ncread(ncdata,'surfvelx'));
	vely=double(ncread(ncdata,'surfvely'));
	vx=InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0);
	vy=InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0);
	vel=sqrt(vx.^2+vy.^2);

	%refine mesh using surface velocities as metric
	md=bamg(md,'hmin',1200,'hmax',15000,'field',vel,'err',5);
	
	name='JKS.Mesh_generation';
	save(name,'md','-v7.3')
end 
if any(steps==2)

	disp('	Step 2: Parameterization');
	md=loadmodel('JKS.Mesh_generation');
	
	md=setmask(md,'','');
	md=parameterize(md,'Jak.par'); 

	%zones of shear margin softening
	weakb=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,'WeakB.exp','node',2);
	pos=find(weakb);
	md.materials.rheology_B(pos)=.3*md.materials.rheology_B(pos);

	name='JKS.Parameterization';
	save(name,'md','-v7.3')
end 
if any(steps==3)

	disp('	Step 3: Control method friction');
	md=loadmodel('JKS.Parameterization');
	md=setflowequation(md,'macayeal','all');

	%Control general
	md.inversion.iscontrol=1;
	md.inversion.nsteps=20;
	md.inversion.step_threshold=0.999*ones(md.inversion.nsteps,1);
	md.inversion.maxiter_per_step=15*ones(md.inversion.nsteps,1);
	md.verbose=verbose('solution',true,'control',true);

	%Cost functions
	md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1);
	md.inversion.cost_functions(1:md.inversion.nsteps,1)=103;
	md.inversion.cost_functions(floor(md.inversion.nsteps/2):md.inversion.nsteps,1)=101;

	%Controls
	md.inversion.control_parameters={'FrictionCoefficient'};
	md.inversion.gradient_scaling(1:md.inversion.nsteps)=50;
	md.inversion.gradient_scaling(5:end)=10;
	md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
	md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);

	%Additional parameters
	md.diagnostic.restol=0.01;
	md.diagnostic.reltol=0.1;
	md.diagnostic.abstol=NaN;
	md.inversion.cost_function_threshold=NaN;

	%Go solve
	md=setcluster(md,cluster);
	md.solver=addoptions(md.solver,NoneAnalysisEnum,asmoptions);
	md.solver=addoptions(md.solver,DiagnosticVertAnalysisEnum,jacobiasmoptions);
	md=solve(md,DiagnosticSolutionEnum);
	
	name='JKS.Control_friction_macayeal';
	save(name,'md','-v7.3')
end 
if any(steps==4)

	disp('	Plotting')
	md=loadmodel('JKS.Control_friction_macayeal');

	plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',...
		'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,...
		'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,...
		'FontSize#all',12,...
		'data',md.initialization.vel,'title','Observed velocity',...
		'data',md.results.DiagnosticSolution.Vel,'title','Modeled Velocity',...
		'colorbar#1','off','colorbar#2','on','colorbartitle#2','[m/yr]',...
		'caxis#1-2',[0,7000],...
		'data',md.geometry.bed,'title','Bed elevation',...
		'data',md.results.DiagnosticSolution.FrictionCoefficient,...
		'title','Friction Coefficient',...
		'colorbar#3','on','colorbartitle#3','[m]','colorbar#4','on');

end

