source: issm/trunk-jpl/src/m/contrib/massbalance/contourmassbalance.m@ 15771

Last change on this file since 15771 was 15771, checked in by Mathieu Morlighem, 12 years ago

CHG: Diagnostic is now Stressbalance

File size: 1.9 KB
Line 
1function dhdt=contourmassbalance(md,file)
2%CONTOURMASSBALANCE - compute the mass balance on a contour
3%
4% Usage:
5% dhdt=contourmassbalance(md,file)
6
7%some checks
8if nargin~=2,
9 help contourmassbalance
10 error('contourmassbalance error message: bad usage');
11end
12if ((length(md.initialization.vx)~=md.mesh.numberofvertices)|(length(md.initialization.vy)~=md.mesh.numberofvertices))
13 error(['thicknessevolution error message: vx and vy should have a length of ' num2str(md.mesh.numberofvertices)])
14end
15if ~exist(file),
16 error(['thicknessevolution error message: file ' file ' not found']);
17end
18
19%Get segments enveloping contour
20segments=contourenvelope(md,file);
21%md.stressbalance.icefront=segments; plotmodel(md,'data','pressureload','expdisp',file);
22
23%get flag list of elements and nodes inside the contour
24nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1);
25elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2));
26
27%conputing Mass flux
28x=md.mesh.x;
29y=md.mesh.y;
30vx=mean(md.initialization.vx(segments(:,1:end-1)),2);
31vy=mean(md.initialization.vy(segments(:,1:end-1)),2);
32H=mean(md.geometry.thickness(segments(:,1:end-1)),2);
33nx=cos(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1)))));
34ny=sin(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1)))));
35L=sqrt((x(segments(:,1))-x(segments(:,2))).^2+(y(segments(:,2))-y(segments(:,1))).^2);
36flux = - md.materials.rho_ice*sum(L.*H.*(vx.*nx+vy.*ny)); %outflux is negative!
37disp(['mass outflux on ' file ' = ' num2str(-flux/10^9) ' Gt/yr']);
38areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
39dhdt=flux/(sum(areas(find(elemin)))*md.materials.rho_ice);
40disp(['dhdt on ' file ' (Flux method) = ' num2str(dhdt) ' m/yr']);
41
42dhdt=thicknessevolution(md);
43in=find(elemin);
44dhdt=sum(dhdt(in).*areas(in))/sum(areas(in));
45disp(['dhdt on ' file ' (divHV method) = ' num2str(dhdt) ' m/yr']);
Note: See TracBrowser for help on using the repository browser.