function md=SetIceSheetBC(md)
%SETICESHEETBC - Create the boundary conditions for diagnostic and thermal models for an IceSheet with no Ice Front
%
%   Usage:
%      md=SetIceSheetBC(md)
%
%   See also: SETICESHELFBC, SETMARINEICESHEETBC

%node on Dirichlet
pos=find(md.mesh.vertexonboundary);
md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
md.diagnostic.spcvx(pos)=0;
md.diagnostic.spcvy(pos)=0;
md.diagnostic.spcvz(pos)=0;
md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);

%Dirichlet Values
if (length(md.inversion.vx_obs)==md.mesh.numberofvertices & length(md.inversion.vy_obs)==md.mesh.numberofvertices)
	disp('      boundary conditions for diagnostic model: spc set as observed velocities');
	md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos);
	md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos);
else
	disp('      boundary conditions for diagnostic model: spc set as zero');
end

%segment on neumann (Ice Front) -> none
if (md.mesh.dimension==2)
	md.diagnostic.icefront=zeros(0,4);
else
	md.diagnostic.icefront=zeros(0,6);
end

%Create zeros basal melting rate and surface mass balance if not specified
if isnan(md.surfaceforcings.precipitation),
	md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
	disp('      no surfaceforcings.precipitation specified: values set as zero');
end
if isnan(md.surfaceforcings.mass_balance),
	md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
	disp('      no surfaceforcings.mass_balance specified: values set as zero');
end
if isnan(md.basalforcings.melting_rate),
	md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
	disp('      no basalforcings.melting_rate specified: values set as zero');
end
if isnan(md.balancethickness.thickening_rate),
	md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
	disp('      no balancethickness.thickening_rate specified: values set as zero');
end

md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);

if (length(md.initialization.temperature)==md.mesh.numberofvertices),
	md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
	pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
	if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
		md.basalforcings.geothermalflux=50*10^-3*ones(md.mesh.numberofvertices,1); %50 mW/m^2
	end
else
	disp('      no thermal boundary conditions created: no observed temperature found');
end
