%Ok, start defining model parameters here

disp('      creating thickness');
md.geometry.surface=2000.-md.mesh.x*tan(0.1*pi/180.); %to have z>0
md.geometry.bed=md.geometry.surface-1000.;
md.geometry.thickness=md.geometry.surface-md.geometry.bed;

disp('      creating drag');
%md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x/2.)).*sin(md.mesh.y*2.*pi/max(md.mesh.x/2.)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)));
md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x)).*sin(md.mesh.y*2.*pi/max(md.mesh.x))));
%Take care of iceshelves: no basal drag
pos=find(md.mask.elementonfloatingice);
md.friction.coefficient(md.mesh.elements(pos,:))=0.;
md.friction.p=ones(md.mesh.numberofelements,1);
md.friction.q=zeros(md.mesh.numberofelements,1);

disp('      creating flow law parameter');
md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);

disp('      boundary conditions for stressbalance model:');
%Create node on boundary first (because we can not use mesh)
md=SetIceSheetBC(md);