[2522] | 1 | function dhdt=contourmassbalance(md,file)
|
---|
[13006] | 2 | %CONTOURMASSBALANCE - compute the mass balance on a contour
|
---|
[2522] | 3 | %
|
---|
| 4 | % Usage:
|
---|
| 5 | % dhdt=contourmassbalance(md,file)
|
---|
| 6 |
|
---|
| 7 | %some checks
|
---|
[2541] | 8 | if nargin~=2,
|
---|
| 9 | help contourmassbalance
|
---|
| 10 | error('contourmassbalance error message: bad usage');
|
---|
| 11 | end
|
---|
[9725] | 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)])
|
---|
[2522] | 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);
|
---|
[9679] | 21 | %md.diagnostic.icefront=segments; plotmodel(md,'data','pressureload','expdisp',file);
|
---|
[2522] | 22 |
|
---|
| 23 | %get flag list of elements and nodes inside the contour
|
---|
[9734] | 24 | nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1);
|
---|
[9733] | 25 | elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2));
|
---|
[2522] | 26 |
|
---|
| 27 | %conputing Mass flux
|
---|
[9734] | 28 | x=md.mesh.x;
|
---|
| 29 | y=md.mesh.y;
|
---|
[9684] | 30 | vx=mean(md.initialization.vx(segments(:,1:end-1)),2);
|
---|
| 31 | vy=mean(md.initialization.vy(segments(:,1:end-1)),2);
|
---|
[9691] | 32 | H=mean(md.geometry.thickness(segments(:,1:end-1)),2);
|
---|
[2522] | 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);
|
---|
[9636] | 36 | flux = - md.materials.rho_ice*sum(L.*H.*(vx.*nx+vy.*ny)); %outflux is negative!
|
---|
[2522] | 37 | disp(['mass outflux on ' file ' = ' num2str(-flux/10^9) ' Gt/yr']);
|
---|
[9734] | 38 | areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
|
---|
[9636] | 39 | dhdt=flux/(sum(areas(find(elemin)))*md.materials.rho_ice);
|
---|
[2522] | 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']);
|
---|