%Start defining model parameters here

%dynamics
md.dt=1*md.yts; %1 year
md.ndt=md.dt*10; 
md.artificial_diffusivity=1;
md.eps_abs=5000;
md.min_thermal_constraints=1;

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

disp('      creating drag');
md.drag_type=2; %0 none 1 plastic 2 viscous
md.drag=200*ones(md.numberofgrids,1); %q=1.
%Take care of iceshelves: no basal drag
pos=find(md.elementoniceshelf);
md.drag(md.elements(pos,:))=0;
md.p=ones(md.numberofelements,1);
md.q=ones(md.numberofelements,1);

disp('      creating temperature');
md.observed_temperature=(273-20)*ones(md.numberofgrids,1);

disp('      creating flow law paramter');
md.B=paterson(md.observed_temperature);
md.n=3*ones(md.numberofelements,1);

%Deal with boundary conditions:
md=SetMarineIceSheetBC(md,'Front.exp');

%Parallel options
md.np=3;
md.time=50;
md.waitonlock=1;