source: issm/oecreview/Archive/21724-22754/ISSM-22152-22153.diff@ 22755

Index: ../trunk-jpl/src/m/geometry/VolumeAboveFloatation.m
===================================================================
--- ../trunk-jpl/src/m/geometry/VolumeAboveFloatation.m (nonexistent)
+++ ../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);