%Start defining model parameters here

% Set up bed topography and ice geometry for a tilted 500m thick slab
    md.geometry.base = .02*md.mesh.x;
    md.geometry.bed = md.geometry.base;
    md.geometry.surface = .02*md.mesh.x + 500;
    md.geometry.thickness = md.geometry.surface - md.geometry.bed;

% Define ice sliding velocity (m/yr) 
    md.initialization.vx = 10^-6*md.constants.yts*ones(md.mesh.numberofvertices,1);
    md.initialization.vy = zeros(md.mesh.numberofvertices,1);
    md.initialization.pressure=zeros(md.mesh.numberofvertices,1);

md.initialization.pressure=zeros(md.mesh.numberofvertices,1);

% Materials
    % Ice flow law parameter (note that the standard parameter A=B^(-3))
    md.materials.rheology_B= (5e-25)^(-1/3)*ones(md.mesh.numberofvertices,1);
    md.initialization.temperature=(273)*ones(md.mesh.numberofvertices,1);
    md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);

%Calving
md.calving.calvingrate=zeros(md.mesh.numberofvertices,1);
%md.calving.spclevelset=NaN(md.mesh.numberofvertices,1);

% Friction law and coefficient
    md.friction=frictionsommers(md.friction);
    md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1);


%Numerical parameters
md.stressbalance.viscosity_overshoot=0.0;
md.masstransport.stabilization=1.;
md.thermal.stabilization=1.;
md.verbose=verbose(0);
md.settings.waitonlock=30;
md.stressbalance.restol=0.05;
md.steadystate.reltol=0.05;
md.stressbalance.reltol=0.05;
md.stressbalance.abstol=NaN;
md.timestepping.time_step=1.;
md.timestepping.final_time=3.;

%GIA: 
md.gia.lithosphere_thickness=100.*ones(md.mesh.numberofvertices,1); % in km
md.gia.mantle_viscosity=1.0*10^21*ones(md.mesh.numberofvertices,1); % in Pa.s
md.materials.lithosphere_shear_modulus=6.7*10^10;                   % in Pa
md.materials.lithosphere_density=3.32;                              % in g/cm^-3
md.materials.mantle_shear_modulus=1.45*10^11;                       % in Pa
md.materials.mantle_density=3.34;                                   % in g/cm^-3

%Boundary conditions:
md=SetIceSheetBC(md);

%Change name so that no test have the same name
A=dbstack;
if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end

