source: issm/trunk/test/Validation/ISMIP/TestD/Square.par@ 3768

%Ok, start defining model parameters here

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

disp('      creating drag');
md.drag_type=2; %0 none 1 plastic 2 viscous
md.drag_coefficient=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
%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 flow law paramter');
md.rheology_B=6.8067*10^7*ones(md.numberofgrids,1);
md.rheology_n=3*ones(md.numberofelements,1);

disp('      boundary conditions for diagnostic model: ');
%Create grid on boundary fist (because wi can not use mesh)
md.gridonboundary=zeros(md.numberofgrids,1);
pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
md.gridonboundary(pos)=1;
md.gridondirichlet_diag=zeros(md.numberofgrids,1);
pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
md.dirichletvalues_diag=zeros(md.numberofgrids,2);
md.segmentonneumann_diag=zeros(0,3);
md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature