%This test is a test from the EISMINT for Ice shelves 
% Vincent Rommelaere 1996

%The goal is to test the prognostic model
%md=bamg(model,'domain','DomainOutline.exp','hmax',4550);
md=bamg(model,'domain','DomainOutline.exp','hmax',3000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'macayeal','all');

%Evolution of the ice shelf
md.ndt=500;
md.dt=1;
md.artificial_diffusivity=0; %Better result with no artificial diffusivity

%0launch transient solution
%md=solve(md,'analysis_type','transient');
%FOR NOW:
md.cluster=oshostname();
md.np=8;
i=0;
time=0;
md.dummy=struct();
pos=find(md.y>199999.9);
connectivity=full(sparse(md.elements(:),1,1));
while(time<500),
	disp(['step ' num2str(i) '/' num2str(500/md.dt)]);
	i=i+1;
	time=time+md.dt;
	md=solve(md,'analysis_type','prognostic2','package','ice');
	%thickness=full(sparse(reshape(md.elements',3*md.numberofelements,1),1,md.results.prognostic2.thickness));
	%md.thickness=thickness./connectivity;
	md.thickness=md.results.prognostic2.thickness;
	md.thickness(pos)=500+100*sin(2*pi*time/200)*ones(size(pos,1),1);
	md.surface=md.bed+md.thickness;
	md.dummy(i).thickness=md.thickness;
	%plotmodel(md,'data',md.thickness,'sectionvalue','CrossLine.exp')
end

%plot results
plotmodel(md,'data',[md.thickness 500+100*sin(2*pi/200*(500-md.y/400))],'sectionvalue','CrossLine.exp')

%Don't forget to add these lines in icetransient2d.m, just before the computation of the thickness to change the thickness on the upper boundary condition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%pos=find(m_p.ys~=0);
%m_p.ys(pos)=500+500/5*sin(2*pi*time/(200*md.yts))*ones(size(pos,1),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%