%Test on the diagnostic model and the prognostic in 2d

%tests 3 and 4: using Glen's flow law
md=model;
md=mesh(md,'DomainOutline.exp',5100); %test3
%md=mesh(md,'DomainOutline.exp',3550); %test4
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'macayeal','all'); %MacAyeal's model and 2d

%Compute solution for MacAyeal's model 
md=solve(md,'diagnostic','ice');

%plot results
plotmodel(md,'data','vx','data','vy','contourlevels',{0,20,40,-20,-40},'contourlevels',{-100,-200,-300,-400,-500,-600})
