source: issm/trunk/test/Par/SquareSheetShelf.par@ 17806

%Start defining model parameters here

%Geometry
hmin=300.;
hmax=1000.;
ymin=min(md.mesh.y);
ymax=max(md.mesh.y);
xmin=min(md.mesh.x);
xmax=max(md.mesh.x);
md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin);
md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
bed_sheet=-md.materials.rho_ice/md.materials.rho_water*(hmax+(hmin-hmax)*(ymax/2-ymin)/(ymax-ymin));
pos=find(md.mesh.y<=ymax/2.);
md.geometry.base(pos)=bed_sheet;
md.geometry.surface=md.geometry.base+md.geometry.thickness;

%Initial velocity 
x     = transpose(ncread('../Data/SquareSheetShelf.nc','x'));
y     = transpose(ncread('../Data/SquareSheetShelf.nc','y'));
vx    = transpose(ncread('../Data/SquareSheetShelf.nc','vx'));
vy    = transpose(ncread('../Data/SquareSheetShelf.nc','vy'));
index = transpose(ncread('../Data/SquareSheetShelf.nc','index'));
md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y);
md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y);
clear vx vy x y index;
md.initialization.vz=zeros(md.mesh.numberofvertices,1);
md.initialization.pressure=zeros(md.mesh.numberofvertices,1);

%Materials
md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1);
md.materials.rheology_B=paterson(md.initialization.temperature);
md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);

%Accumulation and melting
md.surfaceforcings.mass_balance=10.*ones(md.mesh.numberofvertices,1);
md.basalforcings.melting_rate=5.*ones(md.mesh.numberofvertices,1);

%Friction
md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1);
md.friction.coefficient(find(md.mask.groundedice_levelset<0.))=0.;
md.friction.p=ones(md.mesh.numberofelements,1);
md.friction.q=ones(md.mesh.numberofelements,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.;

%Deal with boundary conditions:
md=SetMarineIceSheetBC(md,'../Exp/SquareFront.exp');

%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