function V = VolumeAboveFloatation(md)
%VOLUMEABOVEFLOATATION - returns volume above floatation
%
%   Usage:
%      V = VolumeAboveFloatation(md)

%Special case if 3d
if isa(md.mesh,'mesh3dprisms')
	index = md.mesh.elements2d;
	x = md.mesh.x2d;
	y = md.mesh.y2d;
elseif isa(md.mesh,'mesh2d'),
	index = md.mesh.elements;
	x = md.mesh.x;
	y = md.mesh.y;
else
	error('not supported yet');
end
%1. get some parameters
rho_ice   = md.materials.rho_ice;
rho_water = md.materials.rho_water;

%2. compute averages
base        = mean(md.geometry.base(index),2);
surface     = mean(md.geometry.surface(index),2);
bathymetry  = mean(md.geometry.bed(index),2);

%3. get areas of all triangles
areas = GetAreas(index,x,y);

%4. Compute volume above floatation
V = areas.*(surface-base+min(rho_water/rho_ice*bathymetry,0.));

%5. take out the ones that are outside of levelset or floating
pos = find(min(md.mask.ice_levelset(index),[],2)>0 | min(md.mask.ocean_levelset(index),[],2)<0);
V(pos) = 0;

%sum individual contributions
V = sum(V);
