%Ok, start defining model parameters here disp(' creating thickness'); md.thickness=10*ones(md.numberofnodes,1); md.bed=zeros(md.numberofnodes,1); md.surface=md.bed+md.thickness; disp(' creating drag'); md.friction.coefficient=20*ones(md.numberofnodes,1); %q=1. %no drag is specified in the analytical solution md.friction.p=ones(md.numberofelements,1); md.friction.q=ones(md.numberofelements,1); disp(' creating temperatures'); tmin=238.15; %K st=1.67*10^-2/1000; %k/m; radius=sqrt((md.x).^2+(md.y).^2); md.temperature=(tmin+st*radius); md.basalforcings.geothermalflux=4.2*10^-2*ones(md.numberofnodes,1); disp(' creating flow law paramter'); md.materials.rheology_B=6.81*10^(7)*ones(md.numberofnodes,1); %to have the same B as the analytical solution md.materials.rheology_n=3*ones(md.numberofelements,1); disp(' creating surface mass balance'); smb_max=0.5; %m/yr sb=10^-2/1000; %m/yr/m rel=450*1000; %m md.surfaceforcings.mass_balance=min(smb_max,sb*(rel-radius)); disp(' creating velocities'); constant=0.3; md.vx_obs=constant/2*md.x.*(md.thickness).^-1; md.vy_obs=constant/2*md.y.*(md.thickness).^-1; md.vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2)); md.vx=zeros(md.numberofnodes,1); md.vy=zeros(md.numberofnodes,1); md.vz=zeros(md.numberofnodes,1); md.pressure=zeros(md.numberofnodes,1); %Deal with boundary conditions: disp(' boundary conditions for diagnostic model: '); md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp'); radius=sqrt((md.x).*md.x+(md.y).*md.y); pos=find(radius==min(radius)); md.x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center md.diagnostic.spcvx(pos)=0; md.diagnostic.spcvy(pos)=0; md.diagnostic.spcvz(pos)=0; %parallel options md.timestepping.final_time=50000; %Constants md.materials.rho_ice=910; md.materials.thermalconductivity=2.1; md.materials.latentheat=3.35*10^5; md.materials.beta=8.66*10^-4/(md.materials.rho_ice*md.constants.g); %conversion from K/m to K/Pa md.constants.yts=31556926;