source: issm/trunk/test/Validation/ISMIP/TestD/Square.par@ 16

%Ok, start defining model parameters here

%material parameters
        md.g=9.8;
        md.rho_ice=910;
        md.rho_water=1023;
        di=md.rho_ice/md.rho_water;
        md.yts=365*24*3600;
        md.heatcapacity=2009;
        md.thermalconductivity=2.2; %W/mK
        md.beta=9.8*10^-8;

%Solution parameters
        %parallelization 
        md.cluster='none';
        md.np=2;
        md.time=1;
        md.exclusive=0;

        %statics
        md.lowmem=1;
        md.eps_abs=NaN;
        md.eps_rel=0.01;
        md.penalty_offset=3;
        md.penalty_melting=10^7;
        if md.numberofgrids<1000000,
        md.sparsity=.001;
        else
        md.sparsity=.0001;
        end

        %dynamics
        md.dt=1*md.yts; %1 year
        md.ndt=md.dt*10; 
        md.artificial_diffusivity=1;
        md.viscosity_overshoot=0;

        %control
        md.control_type={'drag'}; %'drag', 'B'
        md.nsteps=5;
        md.tolx=10^-4;
        md.maxiter=20;
        md.optscal=10;
        md.fit='logarithmic'; %'absolute','relative','logarithmic'
        md.meanvel=1000/md.yts; %1000 meters/year
        md.epsvel=eps;


        disp('      creating thickness');
        md.surface=2000-md.x*tan(0.1*pi/180); %to have z>0
        md.bed=md.surface-1000;
        md.thickness=md.surface-md.bed;
        md.firn_layer=0*ones(md.numberofgrids,1);

        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=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));

        %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 temperature');
        md.observed_temperature=(273-20)*ones(md.numberofgrids,1);

        disp('      creating flow law paramter');
        md.B=6.8067*10^7*ones(md.numberofgrids,1);
        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:
        
        %Create grid on boundary fist (because wi can not use mesh)
        md.gridonboundary=zeros(md.numberofgrids,1);
        pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
        md.gridonboundary(pos)=1;

        disp('      boundary conditions for diagnostic model: ');
        %Build gridonicefront, array of boundary grids belonging to the icefront:
        gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
        gridonicefront=double(md.gridonboundary & gridinsideicefront);

        md.gridondirichlet_diag=zeros(md.numberofgrids,1);
        pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
        md.dirichletvalues_diag=zeros(md.numberofgrids,2);

        %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
        %md.segmentonneumann_diag=md.segments(pos,:);
        md.segmentonneumann_diag=zeros(0,3);
        %md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)

        disp('      boundary conditions for prognostic model: ');
        md.gridondirichlet_prog=zeros(md.numberofgrids,1);
        md.dirichletvalues_prog=zeros(md.numberofgrids,1);
        %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
        %md.segmentonneumann_prog=md.segments(pos,:);
        md.segmentonneumann_prog=zeros(0,3);
        md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
        md.neumannvalues_prog(:)=NaN; %free radiation
        
        %pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
        %md.segmentonneumann_prog2=md.segments(pos,:);
        md.segmentonneumann_prog2=zeros(0,3);
        md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
        md.neumannvalues_prog2(:)=NaN; %free radiation
        
        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,:))=50*10^-3; %50 mW/m^2


        

% Some Cielo code, ignore.
if strcmp(md.cluster,'yes')
        ServerDisconnect;
end