source: issm/trunk/src/m/classes/public/modeldefault/defaultparams.m@ 5024

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Some hardcoded parameters for this deck%%%%%%%%%%

%some parameterization for this parameter file :)
thicknesspath=options.thicknesspath;
firnpath=options.firnpath;
surfacepath=options.surfacepath;
mosaicpath=options.mosaicpath;
temperaturepath=options.temperaturepath;
heatfluxpath=options.heatfluxpath;

%Solution parameters

        %parallelization 
        md.cluster=oshostname();
        md.np=8;
        md.time=60;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        disp('      reading bedmap thicknesses');
        md.thickness=InterpFromFile(md.x,md.y,thicknesspath,10);

        disp('      reading firn layer');
        md.firn_layer=InterpFromFile(md.x,md.y,firnpath,0);

        disp('      reading Bamber dem');
        md.surface=InterpFromFile(md.x,md.y,surfacepath,10);

        %correct the surface by taking into account the firn layer
        rho_ice=917;
        rho_firn=830;
        firn_layer_correction=md.firn_layer*(1-rho_firn/rho_ice);
        md.surface=md.surface-firn_layer_correction;
        
        %Correct thickness by assuming hydrostatic equilibrium 10km down the grounding line
        md=ThicknessCorrection(md);

        %some corrections
        minsurf=1;
        pos=find(isnan(md.surface) | (md.surface<=0));
        md.surface(pos)=minsurf;
        pos=find(isnan(md.thickness) | (md.thickness<=0));
        md.thickness(pos)=minsurf/(1-md.rho_ice/md.rho_water);
        md.bed=md.surface-md.thickness;

        disp('      reading velocities from Rignot');
        md=plugvelocities(md,mosaicpath,0);
        
        %drag md.drag or stress
        md.drag_type=2; %0 none 1 plastic 2 viscous
        md.drag_coefficient=300*ones(md.numberofgrids,1); %q=1.
                
        %zones of high md.drag
        %[rhighmd.drag_coefficient]=ArgusContourToMesh(md.elements,md.x,md.y,'HighDrag.exp','node');
        %pos=find(highmd.drag);md.drag_coefficient(pos)=10^3;

        %Take care of iceshelves: no drag md.drag_coefficient
        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);

        %Load md.temperature from Giovinetto:
        disp('      loading temperature');
        md.temperature=InterpFromFile(md.x,md.y,temperaturepath,253);
        while ~isempty(find(isnan(md.temperature))),
                pos=find(isnan(md.temperature));
                if ((pos+1)<=length(md.temperature)),
                        md.temperature(pos)=md.temperature(pos+1);
                else
                        md.temperature(pos)=md.temperature(pos-1);
                end
        end

        %flow law
        disp('      creating flow law paramters');
        md.rheology_n=3*ones(md.numberofelements,1);
        md.rheology_B=paterson(md.temperature);

        %zones of shear margin softening
        %[rweakb]=ArgusContourToMesh(md.elements,md.x,md.y,'Weakmd.BPIG.exp','node');
        %pos=find(weakb);md.rheology_B(pos)=.3*md.rheology_B(pos);
                
        %rifts: none for now.
        %if isstruct(rifts), %we have rifts, choose whether they are full of water or air.
        %       for i=1:length(rifts), 
        %               rifts(i).fill='water'; %choice:'air','water','ice'
        %               rifts(i).friction=10^11;
        %       end
        %end

        disp('      creating accumulation rates');
        md.accumulation_rate=ones(md.numberofgrids,1); %1m/a
        
        %Deal with boundary conditions:

        disp('      thermal model');
        md.melting_rate=zeros(md.numberofgrids,1);
        md.observed_temperature=md.temperature;
        
        disp('      reading geothermal flux');
        load(heatfluxpath);
        md.geothermalflux=InterpFromGridToMesh(x_m,y_m,heatflux_Antarctica,md.x,md.y,80);
        pos=find(md.geothermalflux==0);md.geothermalflux(pos)=80;
        md.geothermalflux=md.geothermalflux/1000; %map is given in mW/m^2, we need it in W/m^2

        disp('      boundary conditions');
        md=SetMarineIceSheetBC(md);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

md.counter=3;