steps=[1:5];

%Hard coded parameters
cluster=generic('name',oshostname,'np',2);
modeldatapath=['../../Data'];

%Run Steps
%Mesh{{{1
if any(steps==1)

	%Mesh parameters
	domain =['./DomainOutline.exp'];
	hmax=400000;
	hmin=80000;
	gradation=1.7;
	err=[15 200];

	%Initial mesh creation
	disp('   Interpolating fields'); 
	md=bamg(model,'domain',domain,'hmax',hmin,'splitcorner',1,'KeepVertices',0);

	ncdata=[modeldatapath '/Greenland_5km_dev1.2.nc'];
	%Use ncdisp('file') to see an ncdump
	x1=double(ncread(ncdata,'x1'));
	y1=double(ncread(ncdata,'y1'));

	velx=double(ncread(ncdata,'surfvelx'));
	vely=double(ncread(ncdata,'surfvely'));

	md.inversion.vx_obs=InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0);
	md.inversion.vy_obs=InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0);
	md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);

	thickness=double(ncread(ncdata,'thk'));
	md.geometry.thickness=InterpFromGridToMesh(x1,y1,thickness',...
		md.mesh.x,md.mesh.y,0);

	disp('   Mesh adaptation'); 
	md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',...
		[md.inversion.vel_obs md.geometry.thickness],'err',err,...
		'splitcorner',1,'KeepVertices',0);
	[md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,+1,39,71);

	name='./Greenland.Mesh_generation';
	save(name,'md');
end% }}}
%Parameterization{{{1
if any(steps==2)

	md=loadmodel('./Greenland.Mesh_generation');

	md=setmask(md,'','');
	md=parameterize(md,'./Greenland.par');

	name='./Greenland.Parameterization'; 
	save(name,'md');
end%}}}
%Control_drag{{{1
if any(steps==3)

	md=loadmodel('./Greenland.Parameterization');
	md=setflowequation(md,'macayeal','all');

	md.diagnostic.restol=0.01;
	md.diagnostic.reltol=0.1;

	md=parametercontroldrag(md,'nsteps',50);
	md.inversion.step_threshold=.99*ones(md.inversion.nsteps,1);
	md.inversion.gradient_scaling(:)=100;

	md.inversion.cost_functions=[];
	md.inversion.cost_functions(1:floor(md.inversion.nsteps),1)=103;
	md.inversion.cost_functions=[103*ones(md.inversion.nsteps,1)  ...
		501*ones(md.inversion.nsteps,1)];
	md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,2);
	md.inversion.cost_functions_coefficients(:,2)=8*10^-15;
	md.inversion.control_parameters={'FrictionCoefficient'};
	md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
	md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);

	md=setcluster(md,cluster);
	md.inversion.iscontrol=1;
	md.verbose=verbose('solution',true,'control',true);
	md=solve(md,DiagnosticSolutionEnum);

	md.friction.coefficient=md.results.DiagnosticSolution.FrictionCoefficient;
	md.initialization.vx=md.results.DiagnosticSolution.Vx;
	md.initialization.vy=md.results.DiagnosticSolution.Vy;

	name='./Greenland.Control_drag'; 
	save(name,'md');
end%}}}
%Transient{{{1
if any(steps==4)

	md=loadmodel('./Greenland.Control_drag');

	md=setflowequation(md,'macayeal','all');
	md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);

	ncdata=[modeldatapath '/Greenland_5km_dev1.2.nc'];
	x1=double(ncread(ncdata,'x1'));
	y1=double(ncread(ncdata,'y1'));

	smb=double(ncread(ncdata,'smb'));
	smb=InterpFromGridToMesh(x1,y1,smb',md.mesh.x,md.mesh.y,0)*...
		1000/md.materials.rho_ice;
	smb=[smb smb smb-1.0];
	md.surfaceforcings.mass_balance=[smb; [1 10 20]];

	md.thermal.spctemperature=[md.initialization.temperature;1]; 

	md.timestepping.time_step=0.2;
	md.timestepping.final_time=20;
	md.settings.output_frequency=1;

	md.inversion.iscontrol=0;
	md.transient.requested_outputs=IceVolumeEnum();

	md=setcluster(md,cluster);
	md.verbose=verbose('solution',true,'module',true,'convergence',true);

	md=solve(md,TransientSolutionEnum);

	name='./Greenland.Transient'; 
	save(name,'md');
end% }}}
%Plot{{{1
if any(steps==5)

	md=loadmodel('./Greenland.Transient');

	%Planview plots
	plotmodel(md,'data',md.results.TransientSolution(end).Vel,'caxis',[1e-1 5000],...
		'log', 10, 'title','Velocity (m/y)', ...
		'data', md.results.TransientSolution(1).SurfaceforcingsMassBalance, ...
		'title', 'Surface mass balance (m/y)', ...
		'data',md.results.TransientSolution(end).Thickness,...
		'title','Thickness (m)', ...
		'data',md.results.TransientSolution(end).Surface, ...
		'title', 'Surface (m)');

	%Line Plots

	figure

	%Plot surface mass balance
	surfmb=[]; for i=1:100; surfmb=[surfmb ...
		md.results.TransientSolution(i).SurfaceforcingsMassBalance]; end
	subplot(3,1,1); plot([0.2:0.2:20],mean(surfmb)); title('Mean Surface mass balance');

	%Plot velocity
	vel=[]; for i=1:100; vel=[vel md.results.TransientSolution(i).Vel]; end
	subplot(3,1,2); plot([0.2:0.2:20],mean(vel)); title('Mean Velocity');

	%Plot Volume
	volume=[]; for i=1:100; volume=[volume md.results.TransientSolution(i).IceVolume]; end
	subplot(3,1,3); plot([0.2:0.2:20],volume); title('Ice Volume');
	xlabel('years')
end% }}}
