% Create model;

update_archive=0;

md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');

%First test macayeal: 
md=setelementstype(md,'macayeal','all');
error;

md.cluster='wilkes';

md=solve(md,'diagnostic_horiz','cielo'); pause(10);
md=loadresultsfromcluster(md,'diagnostic_horiz');
vel_cielo=md.vel;

if update_archive,
	vel_nightly=vel_cielo;
	save diagnostic_horiz_macayeal  vel_nightly
end

load diagnostic_horiz_macayeal 

plotmodel(md,'data',vel_nightly,'title','diagnostic_horiz_macayeal','data',vel_cielo,'data',vel_nightly-vel_cielo,'data',(vel_nightly-vel_cielo)./vel_nightly*100,'figure',1);


%Then test full pattyn

md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=extrude(md,3,3);

md=setelementstype(md,'pattyn','all');

md.cluster='wilkes';

md=solve(md,'diagnostic_horiz','cielo'); pause(10);
md=loadresultsfromcluster(md,'diagnostic_horiz');
vel_cielo=md.vel;

if update_archive,
	vel_nightly=vel_cielo;
	save diagnostic_horiz_pattyn  vel_nightly
end

load diagnostic_horiz_pattyn

plotmodel(md,'data',vel_nightly,'title','diagnostic_horiz_pattyn','data',vel_cielo,'data',vel_nightly-vel_cielo,'data',(vel_nightly-vel_cielo)./vel_nightly*100,'figure',2);

%Then test mixed macayeal pattyn

md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=extrude(md,3,3);

md=setelementstype(md,'pattyn','Pattyn.exp','fill','macayeal');

md.cluster='wilkes';

md=solve(md,'diagnostic_horiz','cielo'); pause(10);
md=loadresultsfromcluster(md,'diagnostic_horiz');
vel_cielo=md.vel;

if update_archive,
	vel_nightly=vel_cielo;
	save diagnostic_horiz_mixed  vel_nightly
end

load diagnostic_horiz_mixed

plotmodel(md,'data',vel_nightly,'title','diagnostic_horiz_mixed','data',vel_cielo,'data',vel_nightly-vel_cielo,'data',(vel_nightly-vel_cielo)./vel_nightly*100,'figure',3,'layer#all',1);

%Now test macayeal with rifts
md=model;
md=mesh(md,'DomainOutlineRifts.exp','Rifts.exp',50000);
md=rifttipsrefine(md,5000,80000);
md=meshprocessrifts(md);

md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'macayeal','all');


md.cluster='wilkes';

md=solve(md,'diagnostic_horiz','cielo'); pause(10);
md=loadresultsfromcluster(md,'diagnostic_horiz');
vel_cielo=md.vel;


if update_archive,
	vel_nightly=vel_cielo;
	save diagnostic_horiz_rifts  vel_nightly
end

load diagnostic_horiz_rifts

s=shearnorm2d(md);
plotmodel(md,'data','vel','title','diagnostic_horiz_rifts','data',log10(s),'data','riftpenetration','data','riftvel','figure',4);
plotmodel(md,'data',vel_nightly,'data',vel_cielo,'data',vel_nightly-vel_cielo,'data',(vel_nightly-vel_cielo)./vel_nightly*100,'figure',5);

%Now test thermal
md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=extrude(md,10,3);
md=setelementstype(md,'pattyn','all');

md.cluster='wilkes';

md.vx=zeros(md.numberofgrids,1);
md.vy=zeros(md.numberofgrids,1);
md.vz=zeros(md.numberofgrids,1);
md.pressure=lithostaticpressure(md.rho_ice,md.g,md.surface,md.thickness);

md=solve(md,'thermalsteady','cielo');pause(10);
md=loadresultsfromcluster(md,'thermalsteady');
temperature_cielo=md.temperature;

if update_archive,
	temperature_nightly=temperature_cielo;
	save thermalsteady  temperature_nightly
end

load thermalsteady

plotmodel(md,'data',temperature_nightly,'title','thermalsteady','data',temperature_cielo,'data',temperature_nightly-temperature_cielo,'data',(temperature_nightly-temperature_cielo)./temperature_nightly*100,'figure',6);



% run control method

update_archive=1;

md=model;
md=mesh(md,'DomainOutline.exp',30000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'macayeal','all');

md.cluster='wilkes';
md.np=8;
md.debug=0;

md.drag(:)=500;
pos=find(md.y<200000);
md.drag(pos)=0;

md=solve(md,'diagnostic_horiz','cielo'); pause(10);
md=loadresultsfromcluster(md,'diagnostic_horiz');
md.vx_obs=md.vx;
md.vy_obs=md.vy;

md=solve(md,'control','cielo');pause(60);
md=loadresultsfromcluster(md,'control');
Jfit=md.cont_J;

if update_archive
	save control Jfit
end
load  control
figure(7),plot(md.cont_J,Jfit);

return;
