Index: /issm/trunk-jpl/src/m/geometry/VolumeAboveFloatation.m
===================================================================
--- /issm/trunk-jpl/src/m/geometry/VolumeAboveFloatation.m	(revision 26912)
+++ /issm/trunk-jpl/src/m/geometry/VolumeAboveFloatation.m	(revision 26913)
@@ -1,7 +1,9 @@
-function V = VolumeAboveFloatation(md)
+function V = VolumeAboveFloatation(md,step,flags)
 %VOLUMEABOVEFLOATATION - returns volume above floatation
 %
 %   Usage:
-%      V = VolumeAboveFloatation(md)
+%      V = VolumeAboveFloatation(md)          % uses model fiels alone
+%      V = VolumeAboveFloatation(md,10)       % Will look at step 10 of transient solution
+%      V = VolumeAboveFloatation(md,10,flags) % Will look at step 10 of transient solution, only flaged elements
 
 %Special case if 3d
@@ -17,4 +19,5 @@
 	error('not supported yet');
 end
+
 %1. get some parameters
 rho_ice   = md.materials.rho_ice;
@@ -22,7 +25,25 @@
 
 %2. compute averages
-base        = mean(md.geometry.base(index),2);
-surface     = mean(md.geometry.surface(index),2);
-bathymetry  = mean(md.geometry.bed(index),2);
+if nargin==1
+	base           = mean(md.geometry.base(index),2);
+	surface        = mean(md.geometry.surface(index),2);
+	bathymetry     = mean(md.geometry.bed(index),2);
+	ice_levelset   = md.mask.ice_levelset;
+	ocean_levelset = md.mask.ocean_levelset;
+else
+	if isprop(md.results.TransientSolution(step),'MaskIceLevelset')
+		ice_levelset   = md.results.TransientSolution(step).MaskIceLevelset;
+	else
+		ice_levelset   = md.mask.ice_levelset;
+	end
+   ocean_levelset = md.results.TransientSolution(step).MaskOceanLevelset;
+   base           = mean(md.results.TransientSolution(step).Base(index),2);
+   surface        = mean(md.results.TransientSolution(step).Surface(index),2);
+	if isprop(md.results.TransientSolution(step),'Bed')
+		bathymetry  = mean(md.results.TransientSolution(step).Bed(index),2);
+	else
+		 bathymetry  = mean(md.geometry.bed(index),2);
+	 end
+end
 
 %3. get areas of all triangles
@@ -33,6 +54,11 @@
 
 %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);
+pos = find(min(ice_levelset(index),[],2)>0 | min(ocean_levelset(index),[],2)<0);
 V(pos) = 0;
+
+%In case we are only looking at one portion of the domain...
+if nargin==3
+	V(find(~flags)) = 0;
+end
 
 %sum individual contributions