1 | function 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
|
---|
8 | if nargin~=2,
|
---|
9 | help contourmassbalance
|
---|
10 | error('contourmassbalance error message: bad usage');
|
---|
11 | end
|
---|
12 | if ((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)])
|
---|
14 | end
|
---|
15 | if ~exist(file),
|
---|
16 | error(['thicknessevolution error message: file ' file ' not found']);
|
---|
17 | end
|
---|
18 |
|
---|
19 | %Get segments enveloping contour
|
---|
20 | segments=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
|
---|
24 | nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1);
|
---|
25 | elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2));
|
---|
26 |
|
---|
27 | %conputing Mass flux
|
---|
28 | x=md.mesh.x;
|
---|
29 | y=md.mesh.y;
|
---|
30 | vx=mean(md.initialization.vx(segments(:,1:end-1)),2);
|
---|
31 | vy=mean(md.initialization.vy(segments(:,1:end-1)),2);
|
---|
32 | H=mean(md.geometry.thickness(segments(:,1:end-1)),2);
|
---|
33 | nx=cos(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1)))));
|
---|
34 | ny=sin(atan2((x(segments(:,1))-x(segments(:,2))) , (y(segments(:,2))-y(segments(:,1)))));
|
---|
35 | L=sqrt((x(segments(:,1))-x(segments(:,2))).^2+(y(segments(:,2))-y(segments(:,1))).^2);
|
---|
36 | flux = - md.materials.rho_ice*sum(L.*H.*(vx.*nx+vy.*ny)); %outflux is negative!
|
---|
37 | disp(['mass outflux on ' file ' = ' num2str(-flux/10^9) ' Gt/yr']);
|
---|
38 | areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
|
---|
39 | dhdt=flux/(sum(areas(find(elemin)))*md.materials.rho_ice);
|
---|
40 | disp(['dhdt on ' file ' (Flux method) = ' num2str(dhdt) ' m/yr']);
|
---|
41 |
|
---|
42 | dhdt=thicknessevolution(md);
|
---|
43 | in=find(elemin);
|
---|
44 | dhdt=sum(dhdt(in).*areas(in))/sum(areas(in));
|
---|
45 | disp(['dhdt on ' file ' (divHV method) = ' num2str(dhdt) ' m/yr']);
|
---|