Index: /issm/trunk/src/m/classes/public/contourmassbalance.m
===================================================================
--- /issm/trunk/src/m/classes/public/contourmassbalance.m	(revision 2522)
+++ /issm/trunk/src/m/classes/public/contourmassbalance.m	(revision 2522)
@@ -0,0 +1,41 @@
+function dhdt=contourmassbalance(md,file)
+%CONTOURMASSBALANCE - compute the mass balance on a given profile
+%
+%   Usage:
+%      dhdt=contourmassbalance(md,file)
+
+%some checks
+if (length(md.vx)~=md.numberofgrids)|(length(md.vy)~=md.numberofgrids)
+	error(['thicknessevolution error message: vx and vy should have a length of ' num2str(md.numberofgrids)])
+end
+if ~exist(file),
+	error(['thicknessevolution error message: file ' file ' not found']);
+end
+
+%Get segments enveloping contour
+segments=contourenvelope(md,file);
+%md.pressureload=segments; plotmodel(md,'data','pressureload','expdisp',file);
+
+%get flag list of elements and nodes inside the contour
+nodein=ContourToMesh(md.elements,md.x,md.y,expread(file,1),'node',1);
+elemin=(sum(nodein(md.elements),2)==size(md.elements,2));
+
+%conputing Mass flux
+x=md.x;
+y=md.y;
+vx=mean(md.vx(segments(:,1:end-1)),2);
+vy=mean(md.vy(segments(:,1:end-1)),2);
+H=mean(md.thickness(segments(:,1:end-1)),2);
+nx=cos(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1)))));
+ny=sin(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1)))));
+L=sqrt((x(segments(:,1))-x(segments(:,2))).^2+(y(segments(:,2))-y(segments(:,1))).^2);
+flux = - md.rho_ice*sum(L.*H.*(vx.*nx+vy.*ny)); %outflux is negative!
+disp(['mass outflux on ' file ' = ' num2str(-flux/10^9) ' Gt/yr']);
+areas=GetAreas(md.elements,md.x,md.y);
+dhdt=flux/(sum(areas(find(elemin)))*md.rho_ice);
+disp(['dhdt on ' file ' (Flux  method) = ' num2str(dhdt) ' m/yr']);
+
+dhdt=thicknessevolution(md);
+in=find(elemin);
+dhdt=sum(dhdt(in).*areas(in))/sum(areas(in));
+disp(['dhdt on ' file ' (divHV method) = ' num2str(dhdt) ' m/yr']);