%Test Name: ISMIP6MeltRateTest_HO
md=triangle(model(),'../Exp/Square.exp',90000.);
md=setmask(md,'../Exp/SquareShelf.exp','');
md=parameterize(md,'../Par/SquareSheetShelf.par');
md.initialization.vx(:)=1.;
md.initialization.vy(:)=1.;
md.geometry.thickness(:)=500-md.mesh.x/10000;
md.geometry.bed =-100-md.mesh.x/1000;
md.geometry.base=-md.geometry.thickness*md.materials.rho_ice/md.materials.rho_water;
md.mask.ocean_levelset=md.geometry.thickness+md.materials.rho_water/md.materials.rho_ice*md.geometry.bed;
pos=find(md.mask.ocean_levelset>=0);
md.geometry.base(pos)=md.geometry.bed(pos);
md.geometry.surface=md.geometry.base+md.geometry.thickness;
md=extrude(md,3,1.1);
md=setflowequation(md,'HO','all');

%Set ISMIP6 melt rate parameters
md.basalforcings = basalforcingsismip6(md.basalforcings);
md.basalforcings.basin_id = zeros(md.mesh.numberofelements,1);
yE = mean(md.mesh.y(md.mesh.elements),2);
pos1 = find(yE>=5e5);   md.basalforcings.basin_id(pos1)=1;
pos2 = find(yE<5e5);     md.basalforcings.basin_id(pos2)=2;
md.basalforcings.num_basins = 2;
md.basalforcings.delta_t   = [0.1 0.2];
md.basalforcings.tf_depths = [0 -1000 -2000];
md.basalforcings.gamma_0 = 14477;
md.basalforcings.islocal = 0;

%Build an artificial tf field (for times 0 and 1, 3 depth layers)
temp1a  = 1.*ones(md.mesh.numberofvertices,1); temp1b = 1.5.*ones(md.mesh.numberofvertices,1); A = [temp1a temp1b ; [0. 1.]];
temp2a  = 2.*ones(md.mesh.numberofvertices,1); temp2b = 2.5.*ones(md.mesh.numberofvertices,1); B = [temp2a temp2b ; [0. 1.]];
temp3a  = 3.*ones(md.mesh.numberofvertices,1); temp3b = 3.5.*ones(md.mesh.numberofvertices,1); C = [temp3a temp3b ; [0. 1.]];
D = cell(1,1,3); D(:,:,1)={A}; D(:,:,2)={B}; D(:,:,3)={C};
md.basalforcings.tf = D;

%Boundary conditions:
md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
md.mask.ice_levelset(find(md.mesh.x==max(md.mesh.x)))=0;

%Model conditions
md.transient.isthermal=0;
md.transient.isstressbalance=1;
md.transient.isgroundingline=1;
md.transient.ismasstransport=1;
md.transient.issmb=1;
md.transient.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate','BasalforcingsIsmip6TfShelf'};
md.groundingline.migration='SubelementMigration';
md.groundingline.friction_interpolation='SubelementFriction1';
md.groundingline.melt_interpolation='SubelementMelt1';
md.timestepping.final_time=1.5;
md.timestepping.time_step=0.5;

md.cluster=generic('name',oshostname(),'np',1);
md=solve(md,'Transient');

field_names     ={'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Pressure1','FloatingiceMeltingrate1','ThermalForcing1',...
	   'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Pressure2','FloatingiceMeltingrate2','ThermalForcing2',...
	   'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Pressure3','FloatingiceMeltingrate3','ThermalForcing3'};
field_tolerances={7e-09,8e-09,8e-09,7e-09,6e-08,7e-08,6e-09,8e-10,7e-08,...
	   7e-09,8e-09,8e-09,7e-09,6e-08,7e-08,6e-09,8e-10,7e-08,...
	   7e-09,8e-09,8e-09,7e-09,6e-08,7e-08,6e-09,8e-10,7e-08};
field_values={...
	   (md.results.TransientSolution(1).Base),...
	   (md.results.TransientSolution(1).Surface),...
	   (md.results.TransientSolution(1).Thickness),...
	   (md.results.TransientSolution(1).MaskOceanLevelset),...
	   (md.results.TransientSolution(1).Vx),...
	   (md.results.TransientSolution(1).Vy),...
	   (md.results.TransientSolution(1).Pressure),...
	   (md.results.TransientSolution(1).BasalforcingsFloatingiceMeltingRate),...
	   (md.results.TransientSolution(1).BasalforcingsIsmip6TfShelf),...
	   (md.results.TransientSolution(2).Base),...
	   (md.results.TransientSolution(2).Surface),...
	   (md.results.TransientSolution(2).Thickness),...
	   (md.results.TransientSolution(2).MaskOceanLevelset),...
	   (md.results.TransientSolution(2).Vx),...
	   (md.results.TransientSolution(2).Vy),...
	   (md.results.TransientSolution(2).Pressure),...
	   (md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),...
	   (md.results.TransientSolution(2).BasalforcingsIsmip6TfShelf),...
	   (md.results.TransientSolution(3).Base),...
	   (md.results.TransientSolution(3).Surface),...
	   (md.results.TransientSolution(3).Thickness),...
	   (md.results.TransientSolution(3).MaskOceanLevelset),...
	   (md.results.TransientSolution(3).Vx),...
	   (md.results.TransientSolution(3).Vy),...
	   (md.results.TransientSolution(3).Pressure),...
	   (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),...
	   (md.results.TransientSolution(3).BasalforcingsIsmip6TfShelf),...
	   };
