% This file can be run to check that the current version of pattyn model is giving 
% coherent results. This test deals with an icesheet with icefront for a 3d model. The geometry 
% is square. Just run this file in Matlab, with a properly setup Ice code. 
% The results of this test will indicate if there is a difference between current computations 
% and archived results.

% Errors  between archived results and the current version will get flagged if they are not within 
% a certain tolerance. The current tolerance is 10^-14. If you have good reasons to believe this 
% tolerance should be lowered (for example, if you are running single precision compilers?), feel 
% free to tweak the tolerance variable.

%Initialize log message for nightly runs.
logstring='Verification/test5_icesheet_icefront_3d: ';

% Create model;

tolerance=10^-13;


md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=extrude(md,10,3);
md=setelementstype(md,'pattyn','all');

%Compute solution with Ice model
md=solve(md,'diagnostic','ice');
vel_pattyn=md.vel;



%Load old velocities and compare with the new model
load archive_pattyn
pos=find(archive_pattyn-vel_pattyn);
error_vel=abs(norm((archive_pattyn(pos)-vel_pattyn(pos))./archive_pattyn(pos),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_vel));

if (error_vel>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR.  Results from Pattyn model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from Pattyn model are comform with the archive version'));
end

%Compute temperature with Ice model
md=solve(md,'thermalsteady','ice');
temp_pattyn=md.temperature;
melting_pattyn=md.melting;

%Load old temperature and compare with the new model
load archive_temperature
post=find(archive_temperature-temp_pattyn);
error_temp=abs(norm((archive_temperature(post)-temp_pattyn(post))./archive_temperature(post),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_temp));

if (error_temp>tolerance);
        logoutput(logstring,sprintf('%s\n','    ERROR. Results from Thermal model differ from the archive version'));
else
        logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from Thermal model are comform with the archive version'));
end

%Load old melting and compare with the new model
load archive_melting
post=find(archive_melting-melting_pattyn);
error_melting=abs(norm((archive_melting(post)-melting_pattyn(post))./archive_melting(post),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_melting));

if (error_melting>tolerance);
        logoutput(logstring,sprintf('%s\n','    ERROR. Results from melting model differ from the archive version'));
else
        logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from melting model are comform with the archive version'));
end

%Compute temperature transient with Ice model
md=solve(md,'thermaltransient','ice');
temptrans_pattyn=md.thermaltransient_results(2).temperature;
meltingtrans_pattyn=md.thermaltransient_results(2).melting;

%Load old temperature transient and compare with the new model
load archive_temperature_transient
post2=find(archive_temperature_transient);
error_temp=abs(norm((archive_temperature_transient(post2)-temptrans_pattyn(post2))./archive_temperature_transient(post2),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_temp));

if (error_temp>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR. Results from thermal transient model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from thermal transient model are comform with the archive version'));
end

%Load old melting transient and compare with the new model
load archive_melting_transient
post2=find(archive_melting_transient);
error_melting_transient=abs(norm((archive_melting_transient(post2)-meltingtrans_pattyn(post2))./archive_melting_transient(post2),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_melting));

if (error_melting_transient>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR. Results from melting transient model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from melting transient model are comform with the archive version'));
end

%Compute prognostic with Ice model
md=solve(md,'prognostic','ice');
thickness=md.new_thickness;

%Load old temperature transient and compare with the new model
load archive_thickness
posth=find(archive_thickness);
error_thickness=abs(norm((archive_thickness(posth)-thickness(posth))./archive_thickness(posth),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_thickness));

if (error_thickness>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR. Results from prognostic model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from prognostic model are comform with the archive version'));
end

