source: issm/trunk/test/Validation/ThermalTests/Melting/Square.par@ 544

%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=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 temperatures');
        md.observed_temperature=273.15*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');
        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