%Ok, start defining model parameters here %parallelization md.cluster='none'; disp(' creating thickness'); h=1000; md.thickness=h*ones(md.numberofgrids,1); md.firn_layer=10*ones(md.numberofgrids,1); md.bed=-1000*ones(md.numberofgrids,1); 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_coefficient=200*ones(md.numberofgrids,1); %q=1. %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 temperatures'); md.observed_temperature=273.15*ones(md.numberofgrids,1); disp(' creating flow law paramter'); md.rheology_B=paterson(md.observed_temperature); md.rheology_n=3*ones(md.numberofelements,1); disp(' creating accumulation rates'); md.accumulation_rate=ones(md.numberofgrids,1)/md.yts; %1m/a md.melting_rate=0*ones(md.numberofgrids,1)/md.yts; %1m/a %Deal with boundary conditions: disp(' boundary conditions for diagnostic model'); md=SetIceShelfBC(md,'Front.exp'); 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,:))=1*10^-3; %1 mW/m^2