heleneseroussi
Thanks!
We also tried now with you recommendation and got a vales closer the output. we still mean between the values as we need lx and ly. I am posting just for completeness of this topic:
`gl = expread('gl1.exp');
sz = size(gl,2);
flux=0;
for i=1:sz,
x =gl(i).x;
y =gl(i).y;
[h,vx,vy]=calc_properties(md,1,x,y);
sec_flux = my_contourmassbalance_line( h, vx, vy,x,y);
flux = flux +sec_flux;
end'
with the functions:
function [H,Vx,Vy] =calc_propertiese(md, step,x,y)
H = InterpFromMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, md.results.TransientSolution(step).Thickness, x, y);
Vx = InterpFromMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, md.results.TransientSolution(step).Vx, x, y);
Vy = InterpFromMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y, md.results.TransientSolution(step).Vy, x, y);
and
function [massFlux] =my_contourmassbalance_line( h, vx, vy,x,y)
lx =x(2:end)-x(1:end-1);
ly =y(2:end)-y(1:end-1);
L=sqrt(lx.^2+ly.^2);
Nx=lx./L;
Ny=ly./L;
Vx=(vx(2:end)+vx(1:end-1))/2;
Vy=(vy(2:end)+vy(1:end-1))/2;
H=(h(2:end)+h(1+end-1))/2;
flux_vec = H.*(Vx.*Nx+Vy.*Ny).*L;
flux = 917*sum(flux_vec);
massFlux = flux/10^12;
disp(['mass outflux = ' num2str(massFlux) ' Gt/yr']);
we got for the entire AIS and all grounding line island:
flux =2.0182e+03
close to the ISSM output. Interestingly for my first set of gl(1) the biggest Island so to say I got 1560.5351 Gt/yr.
Does that mean that the segments function, probably only takes the biggest or first island as the others are not connected?
Because I also also tried just by cutting out the AIS along the grounding line and using the outflow function :
i=1;
md.initialization.vy = md.results.TransientSolution(i).Vy;
md.initialization.vx = md.results.TransientSolution(i).Vx;
md.geometry.thickness = md.results.TransientSolution(i).Thickness;
mds=extract(md,md.mask.ocean_levelset>0);
outflux(mds)
Out flux is 1551.399 Gt/yr
ans =
1.5514e+03
First I was not sure why that did not work, but I guess segements only takes connected points?