
% This file can be run to check that the melting in simple conduction is correctly modeled.
% There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
% at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower
% surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed...
% Just run this file in Matlab, with a properly setup Ice code. 

% Create model;

md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=extrude(md,2,1);
md=setelementstype(md,'Pattyn','all');
thermalboundarycondition;

%analytical results
%melting heat = geothermal flux
%Mb*L*rho=G   => Mb=G/L*rho
melting=md.geothermalflux/(md.rho_ice*md.latentheat)*md.yts;

%modeled  results
md=solve(md,'thermalsteady');

%plot results
relative=abs((md.melting-melting)./melting)*100;
relative(find(md.melting==melting))=0;
plotmodel(md,'data','melting','title','modeled melting','data',melting,'title','analytical melting','data',md.melting-melting,'title','absolute error','data',relative,'title','relative error [%]','layer#all',1)
