Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/MismipGLPosition.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/MismipGLPosition.m (revision 21707) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/MismipGLPosition.m (revision 21707) @@ -0,0 +1,29 @@ +%GLy0 : grounding line position @ y=0km +%GLy40 : grounding line position @ y=40km + +function [GLy0 GLy40] = MismipGLPosition(xGL,yGL,t), + + GLy0 = []; + GLy40 = []; + nsteps = length(t); + + for i=1:nsteps, + + xi = xGL(i,:); + yi = yGL(i,:); + + %find GL position at y=0km + pos = find(min(yi)==yi); + GLy0(i) = max(xi(pos)); %max(xi(pos)) Helene, here + + %find GL position at y=40km + pos = find(yi<45 & yi>35); + x = yi(pos); + v = xi(pos); + xq = [38:0.1:42]; + vq = interp1(x,v,xq,'linear'); + pos = find(xq==40); + GLy40(i) = vq(pos); + end + +end Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/generate_plot.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/generate_plot.m (revision 21707) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/generate_plot.m (revision 21707) @@ -0,0 +1,175 @@ +%L0 and L1 {{{ +md0 = loadmodel('./L0_viscous/Transient_Steadystate_Level_0.mat'); +md1 = loadmodel('./L1/L1_viscous/Transient_steadystate.mat'); +md_L1_R5 = loadmodel('./R5/L1_viscous/Transient_steadystate.mat'); +md_L1_R15 = loadmodel('./R15/L1_viscous/Transient_steadystate.mat'); +md_L1_R30 = loadmodel('./R30/L1_viscous/Transient_steadystate.mat'); +md_L1_R40 = loadmodel('./R40/L1_viscous/Transient_steadystate.mat'); + +[ga0 iv0 ivaf0 GLy400 ne0 t0] = ice_evolution(md0); +[gaL1 ivL1 ivafL1 GLy40L1 neL1 tL1] = ice_evolution(md1); +[ga1 iv1 ivaf1 GLy401 ne1 t1] = ice_evolution(md_L1_R5); +[ga2 iv2 ivaf2 GLy402 ne2 t2] = ice_evolution(md_L1_R15); +[ga3 iv3 ivaf3 GLy403 ne3 t3] = ice_evolution(md_L1_R30); +[ga4 iv4 ivaf4 GLy404 ne4 t4] = ice_evolution(md_L1_R40); +% }}} +%L2 {{{ +md_L2_R5 = loadmodel('./R5/L2_viscous/Transient_steadystate.mat'); +md_L2_R15 = loadmodel('./R15/L2_viscous/Transient_steadystate.mat'); +md_L2_R20 = loadmodel('./R20/L2_viscous/Transient_steadystate.mat'); +md_L2_R30 = loadmodel('./R30/L2_viscous/Transient_steadystate.mat'); + +[ga5 iv5 ivaf5 GLy405 ne5 t5] = ice_evolution(md_L2_R5); +[ga6 iv6 ivaf6 GLy406 ne6 t6] = ice_evolution(md_L2_R15); +[ga7 iv7 ivaf7 GLy407 ne7 t7] = ice_evolution(md_L2_R20); +[ga8 iv8 ivaf8 GLy408 ne8 t8] = ice_evolution(md_L2_R30); +% }}} +%L3 {{{ +md_L3_R30 = loadmodel('./L3_viscous/Transient_steadystate.mat'); +[ga9 iv9 ivaf9 GLy409 ne9 t9] = ice_evolution(md_L3_R30); +% }}} + +%L4 {{{ +md_L4_R30 = loadmodel('./L4_viscous/Transient_steadystate.mat'); +[ga10 iv10 ivaf10 GLy4010 ne10 t10] = ice_evolution(md_L4_R30); +% }}} + +% Scaling {{{ +t0=t0/1000; +tL1=tL1/1000; +t1=t1/1000; +t2=t2/1000; +t3=t3/1000; +t4=t4/1000; +t5=t5/1000; +t6=t6/1000; +t7=t7/1000; +t8=t8/1000; +t9=t9/1000; +t10=t10/1000; + +GLy400=GLy400/1000; +GLy40L1=GLy40L1/1000; +GLy401=GLy401/1000; +GLy402=GLy402/1000; +GLy403=GLy403/1000; +GLy404=GLy404/1000; +GLy405=GLy405/1000; +GLy406=GLy406/1000; +GLy407=GLy407/1000; +GLy408=GLy408/1000; +GLy409=GLy409/1000; +GLy4010=GLy4010/1000; + +ivaf0=ivaf0/10^12; +ivafL1=ivafL1/10^12; +ivaf1=ivaf1/10^12; +ivaf2=ivaf2/10^12; +ivaf3=ivaf3/10^12; +ivaf4=ivaf4/10^12; +ivaf5=ivaf5/10^12; +ivaf6=ivaf6/10^12; +ivaf7=ivaf7/10^12; +ivaf8=ivaf8/10^12; +ivaf9=ivaf9/10^12; +ivaf10=ivaf10/10^12; +% }}} + +figure(1) % {{{ +hold on +plot(t0,GLy400,'*',tL1,GLy40L1,'r+',t1,GLy401,'b',t2,GLy402,'g',t3,GLy403,'m',t4,GLy404,'r','LineWidth',1.5) +axis([0 25 365 455]) +legend('L0','L1','L1-R5','L1-R15','L1-R30','L1-R40','Location','southeast') +title('\fontsize{16} Mismip+ steady state - Rmax analysis') +ylabel('GL @ y=40km (km)') +xlabel('t (kyr)') +hold off +% }}} + +figure(2) % {{{ +hold on +plot(t0,GLy400,'*',t5,GLy405,'b',t6,GLy406,'g',t7,GLy407,'m',t8,GLy408,'r','LineWidth',1.5) +axis([0 25 365 455]) +legend('L0','L2-R5','L2-R15','L2-R20','L2-R30','Location','southeast') +title('\fontsize{16} Mismip+ steady state - Rmax analysis') +ylabel('GL @ y=40km (km)') +xlabel('t (kyr)') +hold off +% }}} + +figure(3) % {{{ +res = [2 1]; +res30 = [2 1 0.5 0.25]; +r5 = [GLy401(end) GLy405(end)]; +r15 = [GLy402(end) GLy406(end)]; +r30 = [GLy403(end) GLy408(end) GLy409(end) GLy4010(end)]; +hold on +plot(log(4),GLy400(end),'^','MarkerSize',7,'MarkerFaceColor','b','MarkerEdgeColor','k'); +plot(log(2),GLy40L1(end),'^','MarkerSize',7,'MarkerFaceColor','k','MarkerEdgeColor','k'); +plot(log(res),r5,'^','MarkerSize',7,'MarkerFaceColor','g','MarkerEdgeColor','k'); +plot(log(res),r15,'^','MarkerSize',7,'MarkerFaceColor','y','MarkerEdgeColor','k'); +plot(log(res30),r30,'^','MarkerSize',7,'MarkerFaceColor','r','MarkerEdgeColor','k'); +plot(log(res(1)),GLy404(end),'^','MarkerSize',7,'MarkerFaceColor','m','MarkerEdgeColor','k'); +plot(log(res(2)),GLy407(end),'^','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k'); +axis([log(0.2) log(5) 437 458]) +legend('L0','L1','R5','R15','R30','L1-R40','L2-R20','Location','northeast') +title('\fontsize{15} Mismip+ steady state: GLpos X Rmax') +ylabel('GL @ y=40km (km)') +xlabel('Resolution (km)') +xtickslabel = [0.25 0.5 1 2 4]; +xticks = log(xtickslabel); +ytickslabel = [440 445 450 455]; +yticks = ytickslabel; +ax = gca; +set(ax,'XTick',xticks); +set(ax,'XTickLabel',xtickslabel); +set(ax,'YTick',yticks); +set(ax,'YTickLabel',ytickslabel); +box on +%set(gca,'xscale','log') +hold off +% }}} + +figure(4) % {{{ +res = [2 1]; +res30 = [2 1 0.5 0.25]; +r5 = [ivaf1(end) ivaf5(end)]; +r15 = [ivaf2(end) ivaf6(end)]; +r30 = [ivaf3(end) ivaf8(end) ivaf9(end) ivaf10(end)]; +hold on +plot(log(4),ivaf0(end),'^','MarkerSize',7,'MarkerFaceColor','b','MarkerEdgeColor','k'); +plot(log(2),ivafL1(end),'^','MarkerSize',7,'MarkerFaceColor','k','MarkerEdgeColor','k'); +plot(log(res),r5,'^','MarkerSize',7,'MarkerFaceColor','g','MarkerEdgeColor','k'); +plot(log(res),r15,'^','MarkerSize',7,'MarkerFaceColor','y','MarkerEdgeColor','k'); +plot(log(res30),r30,'^','MarkerSize',7,'MarkerFaceColor','r','MarkerEdgeColor','k'); +plot(log(res(1)),ivaf4(end),'^','MarkerSize',7,'MarkerFaceColor','m','MarkerEdgeColor','k'); +plot(log(res(2)),ivaf7(end),'^','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k'); +axis([log(0.2) log(5) 34.3 37.2]) +legend('L0','L1','R5','R15','R30','L1-R40','L2-R20','Location','northeast') +title('\fontsize{15} Mismip+ steady state: IVAF X Rmax') +ylabel('Ice volume above floatation (1000 x km3)') +xlabel('Resolution (km)') +xtickslabel = [0.25 0.5 1 2 4]; +xticks = log(xtickslabel); +ytickslabel = [34.5 35.0 35.5 36.0 36.5 37.0]; +yticks = ytickslabel; +ax = gca; +set(ax,'XTick',xticks); +set(ax,'XTickLabel',xtickslabel); +set(ax,'YTick',yticks); +set(ax,'YTickLabel',ytickslabel); +box on +%set(gca,'xscale','log') +hold off +% }}} + +figure(5) % {{{ +hold on +plot(t0,GLy400,'*',t3,GLy403,'b',t8,GLy408,'m',t9,GLy409,'r',t10,GLy4010,'g','LineWidth',1.5) +axis([0 25 365 460]) +legend('L0','L1-R30','L2-R30','L3-R30','L4-R30','Location','southeast') +title('\fontsize{16} Mismip+ steady state - Rmax analysis') +ylabel('GL @ y=40km (km)') +xlabel('t (kyr)') +hold off +% }}} Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/gl_position.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/gl_position.m (revision 21707) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/gl_position.m (revision 21707) @@ -0,0 +1,187 @@ +function [glx gly] = gl_position(md,step,level), + + %initialization of some variables + data = md.results.TransientSolution(step).MaskGroundediceLevelset; + index = md.results.TransientSolution(step).MeshElements; + x = md.results.TransientSolution(step).MeshX; + y = md.results.TransientSolution(step).MeshY; + numberofelements = size(index,1); + elementslist = 1:numberofelements; + c = []; + h = []; + + %get unique edges in mesh + %1: list of edges + edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])]; + %2: find unique edges + [edges,I,J]=unique(sort(edges,2),'rows'); + %3: unique edge numbers + vec=J; + %4: unique edges numbers in each triangle (2 triangles sharing the same edge will have + % the same edge number) + edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)]; + + %segments [nodes1 nodes2] + Seg1=index(:,[1 2]); + Seg2=index(:,[2 3]); + Seg3=index(:,[3 1]); + + %segment numbers [1;4;6;...] + Seg1_num=edges_tria(:,1); + Seg2_num=edges_tria(:,2); + Seg3_num=edges_tria(:,3); + + %value of data on each tips of the segments + Data1=data(Seg1); + Data2=data(Seg2); + Data3=data(Seg3); + + %get the ranges for each segment + Range1=sort(Data1,2); + Range2=sort(Data2,2); + Range3=sort(Data3,2); + + %find the segments that contain this value + pos1=(Range1(:,1)level); + pos2=(Range2(:,1)level); + pos3=(Range3(:,1)level); + + %get elements + poselem12=(pos1 & pos2); + poselem13=(pos1 & pos3); + poselem23=(pos2 & pos3); + poselem=find(poselem12 | poselem13 | poselem23); + numelems=length(poselem); + + %if no element has been flagged, skip to the next level + if numelems==0, + return, + end + + %go through the elements and build the coordinates for each segment (1 by element) + x1=zeros(numelems,1); + x2=zeros(numelems,1); + y1=zeros(numelems,1); + y2=zeros(numelems,1); + edge_l=zeros(numelems,2); + + for j=1:numelems, + + weight1=(level-Data1(poselem(j),1))/(Data1(poselem(j),2)-Data1(poselem(j),1)); + weight2=(level-Data2(poselem(j),1))/(Data2(poselem(j),2)-Data2(poselem(j),1)); + weight3=(level-Data3(poselem(j),1))/(Data3(poselem(j),2)-Data3(poselem(j),1)); + + if poselem12(poselem(j)); + + x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1))); + x2(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1))); + y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1))); + y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1))); + edge_l(j,1)=Seg1_num(poselem(j)); + edge_l(j,2)=Seg2_num(poselem(j)); + + elseif poselem13(poselem(j)), + + x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1))); + x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1))); + y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1))); + y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1))); + edge_l(j,1)=Seg1_num(poselem(j)); + edge_l(j,2)=Seg3_num(poselem(j)); + + elseif poselem23(poselem(j)), + + x1(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1))); + x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1))); + y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1))); + y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1))); + edge_l(j,1)=Seg2_num(poselem(j)); + edge_l(j,2)=Seg3_num(poselem(j)); + else + %it shoud not go here + end + end + + %now that we have the segments, we must try to connect them... + + %loop over the subcontours + indice=0; + while ~isempty(edge_l), + indice=indice+1; + + %take the right edge of the second segment and connect it to the next segments if any + e1=edge_l(1,1); e2=edge_l(1,2); + xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)]; + + %erase the lines corresponding to this edge + edge_l(1,:)=[]; + x1(1)=[]; x2(1)=[]; + y1(1)=[]; y2(1)=[]; + + [ro1,co1]=find(edge_l==e1); + + while ~isempty(ro1) + + if co1==1, + xc=[x2(ro1);xc]; yc=[y2(ro1);yc]; + + %next edge: + e1=edge_l(ro1,2); + + else + xc=[x1(ro1);xc]; yc=[y1(ro1);yc]; + + %next edge: + e1=edge_l(ro1,1); + end + + %erase the lines of this + edge_l(ro1,:)=[]; + x1(ro1)=[]; x2(ro1)=[]; + y1(ro1)=[]; y2(ro1)=[]; + + %next connection + [ro1,co1]=find(edge_l==e1); + end + + %same thing the other way (to the right) + [ro2,co2]=find(edge_l==e2); + + while ~isempty(ro2) + + if co2==1, + xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)]; + + %next edge: + e2=edge_l(ro2,2); + else + xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)]; + + %next edge: + e2=edge_l(ro2,1); + end + + %erase the lines of this + edge_l(ro2,:)=[]; + x1(ro2)=[]; x2(ro2)=[]; + y1(ro2)=[]; y2(ro2)=[]; + + %next connection + [ro2,co2]=find(edge_l==e2); + end + + % Update the CS data structure as per "contours.m" + % so that clabel works + c = horzcat(c,[level, xc'; length(xc), yc']); + % y0=find(c(2,:)==0); + % y50=find(c(2,:)==50000); + % gl0(glstep)=c(1,y0); + % gl50(glstep)=c(1,y50); + glx=c(1,1:end)'; + gly=c(2,1:end)'; + pos=find(glx~=0); + glx=glx(pos); + gly=gly(pos); + +end + %min(c(1,2:end)) Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m (revision 21707) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m (revision 21707) @@ -0,0 +1,34 @@ +%ga: grounded area +%iv: ice volume +%ivaf: ice volume above floatation +%GLy40 : grounding line position @ y=40km +%nelem : number of elements + +function [ga iv ivaf GLy40 nelem t] = ice_evolution(md), + + ga = []; + iv = []; + ivaf = []; + GLy40 = []; + nelem = []; + t = []; + nsteps = length(md.results.TransientSolution); + + for i=1:nsteps, + ga(i) = md.results.TransientSolution(i).GroundedArea; + iv(i) = md.results.TransientSolution(i).IceVolume; + ivaf(i) = md.results.TransientSolution(i).IceVolumeAboveFloatation; + nelem(i) = size(md.results.TransientSolution(i).MeshElements,1); + t(i) = md.results.TransientSolution(i).time; + %find GL position at y=40km + [glx gly] = gl_position(md,i,0); + pos=find(gly<45000 & gly > 35000); + x=gly(pos); + v=glx(pos); + xq=[38000:100:42000]; + vq = interp1(x,v,xq,'linear'); + pos=find(xq==40000); + GLy40(i)=vq(pos); + end + +end Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/writeNetCDF.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/writeNetCDF.m (revision 21707) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/writeNetCDF.m (revision 21707) @@ -0,0 +1,112 @@ +function writeNetCDF(md0,md,step,dt,ncfile), + + xGL={}; + yGL={}; + iceThicknessGL={}; + uBaseGL={}; + vBaseGL={}; + iceVolume=[]; + iceVAF=[]; + groundedArea=[]; + time=[]; + + %Inserting time 0. md0 must be last experiment (e.g., Ice1r for Ice1ra) + x = md0.results.TransientSolution(end).MeshX; + y = md0.results.TransientSolution(end).MeshY; + time(1) = 0; + [xgl_step ygl_step] = gl_position(md0,length(md0.results.TransientSolution),0); + xGL{1} = xgl_step; + yGL{1} = ygl_step; + iceVolume(1) = md0.results.TransientSolution(end).IceVolume; + iceVAF(1) = md0.results.TransientSolution(end).IceVolumeAboveFloatation; + groundedArea(1) = md0.results.TransientSolution(end).GroundedArea; + iceThicknessGL{1} = griddata(x,y,md0.results.TransientSolution(end).Thickness,xgl_step,ygl_step); + uBaseGL{1} = griddata(x,y,md0.results.TransientSolution(end).Vx,xgl_step,ygl_step); + vBaseGL{1} = griddata(x,y,md0.results.TransientSolution(end).Vy,xgl_step,ygl_step); + + for i=2:length(step), + x = md.results.TransientSolution(step(i)).MeshX; + y = md.results.TransientSolution(step(i)).MeshY; + time(i)=md.results.TransientSolution(step(i)).time; + [xgl_step ygl_step]=gl_position(md,step(i),0); + xGL{i}=xgl_step; + yGL{i}=ygl_step; + iceVolume(i)=md.results.TransientSolution(step(i)).IceVolume; + iceVAF(i)=md.results.TransientSolution(step(i)).IceVolumeAboveFloatation; + groundedArea(i)=md.results.TransientSolution(step(i)).GroundedArea; + iceThicknessGL{i}=griddata(x,y,md.results.TransientSolution(step(i)).Thickness,xgl_step,ygl_step); + uBaseGL{i}=griddata(x,y,md.results.TransientSolution(step(i)).Vx,xgl_step,ygl_step); + vBaseGL{i}=griddata(x,y,md.results.TransientSolution(step(i)).Vy,xgl_step,ygl_step); + end + uSurfaceGL=uBaseGL; + vSurfaceGL=vBaseGL; + uMeanGL=uBaseGL; + vMeanGL=vBaseGL; + + %Create netcdf + mode = netcdf.getConstant('NETCDF4'); + mode = bitor(mode,netcdf.getConstant('CLASSIC_MODEL')); + ncid=netcdf.create(ncfile,mode); + + %General attributes + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author','Helene Seroussi (helene.seroussi@jpl.nasa.gov)'); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Model','ISSM (Ice Sheet System Model)'); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Notes','Experiments performed at Caltech Jet Propulsion Laboratory, Pasadena'); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',date()); + + %Define dimensions + gl_id = netcdf.defDim(ncid,'nPointGL',netcdf.getConstant('NC_UNLIMITED')); + ntime_id = netcdf.defDim(ncid,'nTime',length(time)); + + %Define variables + volume_var_id = netcdf.defVar(ncid,'iceVolume','NC_FLOAT',[ntime_id]); + VAF_var_id = netcdf.defVar(ncid,'iceVAF','NC_FLOAT',[ntime_id]); + grounded_var_id = netcdf.defVar(ncid,'groundedArea','NC_FLOAT',[ntime_id]); + xgl_var_id = netcdf.defVar(ncid,'xGL','NC_FLOAT',[ntime_id,gl_id]); + ygl_var_id = netcdf.defVar(ncid,'yGL','NC_FLOAT',[ntime_id,gl_id]); + thickness_var_id = netcdf.defVar(ncid,'iceThicknessGL','NC_FLOAT',[ntime_id,gl_id]); + ubase_var_id = netcdf.defVar(ncid,'uBaseGL','NC_FLOAT',[ntime_id,gl_id]); + vbase_var_id = netcdf.defVar(ncid,'vBaseGL','NC_FLOAT',[ntime_id,gl_id]); + usurface_var_id = netcdf.defVar(ncid,'uSurfaceGL','NC_FLOAT',[ntime_id,gl_id]); + vsurface_var_id = netcdf.defVar(ncid,'vSurfaceGL','NC_FLOAT',[ntime_id,gl_id]); + umean_var_id = netcdf.defVar(ncid,'uMeanGL','NC_FLOAT',[ntime_id,gl_id]); + vmean_var_id = netcdf.defVar(ncid,'vMeanGL','NC_FLOAT',[ntime_id,gl_id]); + time_var_id = netcdf.defVar(ncid,'time','NC_FLOAT',[ntime_id]); + + %Define default fill values + if(false), + netcdf.defVarFill(ncid,xgl_var_id,false,NaN); + netcdf.defVarFill(ncid,ygl_var_id,false,NaN); + netcdf.defVarFill(ncid,thickness_var_id,false,NaN); + netcdf.defVarFill(ncid,ubase_var_id,false,NaN); + netcdf.defVarFill(ncid,vbase_var_id,false,NaN); + netcdf.defVarFill(ncid,usurface_var_id,false,NaN); + netcdf.defVarFill(ncid,vsurface_var_id,false,NaN); + netcdf.defVarFill(ncid,umean_var_id,false,NaN); + netcdf.defVarFill(ncid,vmean_var_id,false,NaN); + end + + netcdf.endDef(ncid); + + %Write variables + netcdf.putVar(ncid,volume_var_id,iceVolume); + netcdf.putVar(ncid,VAF_var_id,iceVAF); + netcdf.putVar(ncid,grounded_var_id,groundedArea); + for i=1:length(time), + netcdf.putVar(ncid,xgl_var_id,[i-1,0],[1,length(xGL{i})],xGL{i}'); + netcdf.putVar(ncid,ygl_var_id,[i-1,0],[1,length(xGL{i})],yGL{i}'); + netcdf.putVar(ncid,thickness_var_id,[i-1,0],[1,length(xGL{i})],iceThicknessGL{i}'); + netcdf.putVar(ncid,ubase_var_id,[i-1,0],[1,length(xGL{i})],uBaseGL{i}'); + netcdf.putVar(ncid,vbase_var_id,[i-1,0],[1,length(xGL{i})],vBaseGL{i}'); + netcdf.putVar(ncid,usurface_var_id,[i-1,0],[1,length(xGL{i})],uSurfaceGL{i}'); + netcdf.putVar(ncid,vsurface_var_id,[i-1,0],[1,length(xGL{i})],vSurfaceGL{i}'); + netcdf.putVar(ncid,umean_var_id,[i-1,0],[1,length(xGL{i})],uMeanGL{i}'); + netcdf.putVar(ncid,vmean_var_id,[i-1,0],[1,length(xGL{i})],vMeanGL{i}'); + end + netcdf.putVar(ncid,time_var_id,time+dt); + + %Close netcdf + netcdf.close(ncid) + +end +