Index: /issm/trunk-jpl/src/m/geometry/VolumeAboveFloatation.m
===================================================================
--- /issm/trunk-jpl/src/m/geometry/VolumeAboveFloatation.m	(revision 22153)
+++ /issm/trunk-jpl/src/m/geometry/VolumeAboveFloatation.m	(revision 22153)
@@ -0,0 +1,39 @@
+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.groundedice_levelset(index),[],2)<0);
+V(pos) = 0;
+
+%sum individual contributions
+V = sum(V);