%Start defining model parameters here

disp('      creating thickness');
hmin=300;
hmax=1000;
ymin=min(md.mesh.y);
ymax=max(md.mesh.y);
md.geometry.thickness = hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin);
md.geometry.base      = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
md.geometry.surface   = md.geometry.base+md.geometry.thickness;

disp('      creating drag');
md.friction.coefficient=200*ones(md.mesh.numberofvertices,1);
md.friction.coefficient(find(md.mask.ocean_levelset<0.))=0.;
md.friction.p=ones(md.mesh.numberofelements,1);
md.friction.q=ones(md.mesh.numberofelements,1);

disp('      initial velocity');
md.initialization.vx=zeros(md.mesh.numberofvertices,1);
md.initialization.vy=zeros(md.mesh.numberofvertices,1);
md.initialization.vz=zeros(md.mesh.numberofvertices,1);
md.initialization.vel=zeros(md.mesh.numberofvertices,1);

disp('      creating flow law parameter');
md.materials.rheology_B=paterson((273-20)*ones(md.mesh.numberofvertices,1));
md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
md.damage.D=zeros(md.mesh.numberofvertices,1);

disp('      creating boundary conditions');
md=SetIceShelfBC(md,'Front.exp');
