%Ok, start defining model parameters here

%material parameters
	md.g=9.8;
	md.rho_ice=917;
	md.rho_water=1023;
	di=md.rho_ice/md.rho_water;
	md.yts=365*24*3600;
	md.heatcapacity=2009;
	md.thermalconductivity=2.2; %W/mK
	md.beta=9.8*10^-8;

%Solution parameters
	%parallelization 
	md.cluster='none';
	md.np=2;
	md.time=1;
	md.exclusive=0;

	%statics
	md.eps_rel=0.01;
	md.eps_abs=10;
	md.lowmem=1;
	if md.numberofgrids<1000000,
	md.sparsity=.001;
	else
	md.sparsity=.0001;
	end

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

	%control
	md.control_type={'drag'}; %'drag', 'B'
	md.nsteps=5;
	md.tolx=10^-4;
	md.maxiter=20;
	md.optscal=10;
	md.fit='logarithmic'; %'absolute','relative','logarithmic'
	md.meanvel=1000/md.yts; %1000 meters/year
	md.epsvel=eps;


	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.firn_layer=10*ones(md.numberofgrids,1);
	md.bed=-di*md.thickness;
	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=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);

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

	%Deal with boundary conditions:
	
	disp('      boundary conditions for diagnostic model: ');
	%Build gridonicefront, array of boundary grids belonging to the icefront:
	gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
	gridonicefront=double(md.gridonboundary & gridinsideicefront);

	md.gridondirichlet_diag=zeros(md.numberofgrids,1);
	pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
	md.dirichletvalues_diag=zeros(md.numberofgrids,2);

	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
	md.segmentonneumann_diag=md.segments(pos,:);
	md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)

	disp('      boundary conditions for prognostic model: ');
	md.gridondirichlet_prog=zeros(md.numberofgrids,1);
	md.dirichletvalues_prog=zeros(md.numberofgrids,1);
	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
	md.segmentonneumann_prog=md.segments(pos,:);
	md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
	md.neumannvalues_prog(:)=NaN; %free radiation
	
	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
	md.segmentonneumann_prog2=md.segments(pos,:);
	md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
	md.neumannvalues_prog2(:)=NaN; %free radiation
	
	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,:))=50*10^-3; %50 mW/m^2


	

% Some Cielo code, ignore.
if strcmp(md.cluster,'yes')
	ServerDisconnect;
end