%Ok, start defining model parameters here

%parallelization 
md.cluster='none';

disp('      creating thickness');
h=1000;
md.thickness=h*ones(md.numberofgrids,1);
md.firn_layer=10*ones(md.numberofgrids,1);
md.bed=-1000*ones(md.numberofgrids,1);
md.surface=md.bed+md.thickness;

disp('      creating velocities');
md.vx_obs=zeros(md.numberofgrids,1);
md.vy_obs=zeros(md.numberofgrids,1);
md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);

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

disp('      creating temperatures');
md.observed_temperature=273.15*ones(md.numberofgrids,1);

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

disp('      creating accumulation rates');
md.accumulation_rate=ones(md.numberofgrids,1)/md.yts; %1m/a
md.melting_rate=0*ones(md.numberofgrids,1)/md.yts; %1m/a

%Deal with boundary conditions:

disp('      boundary conditions for diagnostic model');
md=SetIceShelfBC(md,'Front.exp');

disp('      boundary conditions for thermal model');
md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
md.dirichletvalues_thermal=md.observed_temperature;
md.geothermalflux=zeros(md.numberofgrids,1); 
pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2