%Test Name: SquareShelfTranIspddIsdeltaSSA3d
md=triangle(model(),'../Exp/Square.exp',600000.);%180000
md=setmask(md,'all','');
md=parameterize(md,'../Par/SquareShelf.par');

%md.verbose=verbose('all');

% Use of ispdd methods
md.smb = SMBpdd();
md.smb.isdelta18o=0;
md.smb.ismungsm=1;

% time steps and resolution
md.timestepping.time_step=20;
md.settings.output_frequency=1;
md.timestepping.final_time=60;

% creating Present day and lgm temperatures
% Same temperature over the all region:
curve=sin(0:pi/11:pi);
tmonth(1:12)=238.15+20.0+20.0*curve;
for imonth=0:11
    md.smb.temperatures_presentday(1:md.mesh.numberofvertices,imonth+1)=tmonth(imonth+1);
    md.smb.temperatures_lgm(1:md.mesh.numberofvertices,imonth+1)=tmonth(imonth+1)-20.;
end

% creating initialization and spc temperatures initialization and spc
md.thermal.spctemperature=mean(md.smb.temperatures_lgm(1:md.mesh.numberofvertices,1:12),2); %-10*ones(md.mesh.numberofvertices,1);
md.thermal.spctemperature=repmat(md.thermal.spctemperature,1,md.timestepping.final_time/md.timestepping.time_step);
itemp=0:md.timestepping.time_step:md.timestepping.final_time-md.timestepping.time_step;
md.thermal.spctemperature(md.mesh.numberofvertices+1,:)=itemp;

md.initialization.temperature=md.smb.temperatures_lgm(1:md.mesh.numberofvertices,1); %*ones(md.mesh.numberofvertices,1);
md.smb = initialize(md.smb,md);

% creating precipitation
for imonth=0:11
    md.smb.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
    md.smb.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
end

for iint=1:(md.timestepping.final_time/md.timestepping.time_step)+1
    md.smb.Pfac(1,iint)=0.15*iint;
    md.smb.Tdiff(1,iint)=0.15*iint;
    md.smb.sealev(1,iint)=0.15*iint;
    % Year of each data point
    md.smb.Pfac(2,iint)=(iint-1)*20;
    md.smb.Tdiff(2,iint)=(iint-1)*20;
    md.smb.sealev(2,iint)=(iint-1)*20;
end

md.smb.issetpddfac=1;
md.smb.pddfac_snow=2;
md.smb.pddfac_ice=2;

md=extrude(md,3,1);

md=setflowequation(md,'SSA','all');
md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster
md=solve(md,'Transient');

%Fields and tolerances to track changes
field_names     ={'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1','SmbMassBalance1',...
						'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Temperature2','BasalforcingsGroundediceMeltingRate2','SmbMassBalance2',...
						'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','BasalforcingsGroundediceMeltingRate3','SmbMassBalance3'};
field_tolerances={...
	1e-13,1e-13,1e-13,1e-13,1e-13,1e-8,1e-8,1e-8,1e-13,1e-8,1e-13,...
	1e-13,1e-13,1e-13,1e-13,1e-13,1e-8,1e-8,1e-8,1e-13,1e-7,1e-13,...
	1e-13,1e-13,1e-08,1e-13,1e-13,1e-8,1e-8,1e-8,1e-13,5e-7,1e-13};
field_values={...
	(md.results.TransientSolution(1).Vx),...
	(md.results.TransientSolution(1).Vy),...
	(md.results.TransientSolution(1).Vz),...
	(md.results.TransientSolution(1).Vel),...
	(md.results.TransientSolution(1).Pressure),...
	(md.results.TransientSolution(1).Base),...
	(md.results.TransientSolution(1).Surface),...
	(md.results.TransientSolution(1).Thickness),...
	(md.results.TransientSolution(1).Temperature),...
	(md.results.TransientSolution(1).BasalforcingsGroundediceMeltingRate),...
	(md.results.TransientSolution(1).SmbMassBalance),...
	(md.results.TransientSolution(2).Vx),...
	(md.results.TransientSolution(2).Vy),...
	(md.results.TransientSolution(2).Vz),...
	(md.results.TransientSolution(2).Vel),...
	(md.results.TransientSolution(2).Pressure),...
	(md.results.TransientSolution(2).Base),...
	(md.results.TransientSolution(2).Surface),...
	(md.results.TransientSolution(2).Thickness),...
	(md.results.TransientSolution(2).Temperature),...
	(md.results.TransientSolution(2).BasalforcingsGroundediceMeltingRate),...
	(md.results.TransientSolution(2).SmbMassBalance),...
	(md.results.TransientSolution(3).Vx),...
	(md.results.TransientSolution(3).Vy),...
	(md.results.TransientSolution(3).Vz),...
	(md.results.TransientSolution(3).Vel),...
	(md.results.TransientSolution(3).Pressure),...
	(md.results.TransientSolution(3).Base),...
	(md.results.TransientSolution(3).Surface),...
	(md.results.TransientSolution(3).Thickness),...
	(md.results.TransientSolution(3).Temperature),...
	(md.results.TransientSolution(3).BasalforcingsGroundediceMeltingRate),...
	(md.results.TransientSolution(3).SmbMassBalance),...
	};