md=mesh(model,'../Exp/Square.exp',150000);
md=geography(md,'all','');
md=parameterize(md,'../Par/SquareShelf.par');
md=setelementstype(md,'macayeal','all');
md.cluster=none;

%redo the parameter file for this special shelf. 
%constant thickness, constrained (vy=0) flow into an icefront, 
%from 0 m/yr at the grounding line.

%tighten
md.eps_res=10^-4;

%needed later
ymin=min(md.y);
ymax=max(md.y);
xmin=min(md.x);
xmax=max(md.x);

di=md.rho_ice/md.rho_water;

h=1000;
md.thickness=h*ones(md.numberofgrids,1);
md.bed=-md.rho_ice/md.rho_water*md.thickness;
md.surface=md.bed+md.thickness;

%Initial velocity and pressure
md.vx=zeros(md.numberofgrids,1);
md.vy=zeros(md.numberofgrids,1);
md.vz=zeros(md.numberofgrids,1);
md.pressure=zeros(md.numberofgrids,1);

%Materials
md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
md.rheology_B=paterson(md.observed_temperature);
md.rheology_n=3*ones(md.numberofelements,1);
md.temperature=md.observed_temperature;

%Boundary conditions:
md.spcvelocity=zeros(md.numberofgrids,6);

%constrain flanks to 0 normal velocity
pos=find(md.x==xmin | md.x==xmax);
md.spcvelocity(pos,1)=1;
md.spcvelocity(pos,3)=0;

%constrain grounding line to 0 velocity
pos=find(md.y==ymin);
md.spcvelocity(pos,1:2)=1;
md.spcvelocity(pos,3:4)=0;

%icefront
gridonicefront=zeros(md.numberofgrids,1);
pos=find(md.y==ymax); gridonicefront(pos)=1;
pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2))); pressureload=md.segments(pos,:);
pressureload=[pressureload WaterEnum*md.elementoniceshelf(pressureload(:,end))];
md.pressureload=pressureload;

md=solve(md,'analysis_type',DiagnosticSolutionEnum);

%create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292.
%ey_c=(md.rho_ice*md.g*(1-di)*md.thickness./(4*md.rheology_B)).^3;
%vy_c=ey_c.*md.y*md.yts;

%Fields and tolerances to track changes
field_names     ={'Vy'};
field_tolerances={1e-13};
field_values={PatchToVec(md.results.DiagnosticSolution.Vy)};