Index: /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/interpBedmap2.m =================================================================== --- /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/interpBedmap2.m (revision 27069) +++ /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/interpBedmap2.m (revision 27069) @@ -0,0 +1,47 @@ +function [dataout] = interpBedmap2(X,Y,string), +%INTERPBEDMAP2 - interpolate bedmap2 data +% +% Available data: +% 1. bed is bed height +% 2. surface is surface height +% 3. thickness is ice thickness +% 4. icemask_grounded_and_shelves is a mask file showing the grounding line and the extent of the floating ice shelves +% 5. rockmask is a mask file showing rock outcrops +% 6. lakemask_vostok is a mask file showing the extent of the lake cavity of Lake Vostok +% 7. bed_uncertainty is the bed uncertainty grid shown in figure 12 of the manuscript +% 8. thickness_uncertainty_5km is the thickness uncertainty grid shown in figure 11 of the manuscript +% 9. data_coverage is a binary grid showing the dis tribution of ice thickness data used in the grid of ice thickness +% +% Usage: +% [dataout] = interpBedmap2(X,Y,string) + +%reqad data +%path = '/u/astrid-r1b/ModelData/BedMap2/bedmap2_bin/'; +%path = '~/issm-jpl/proj-group/ModelData/BedMap2/bedmap2_bin/'; +path = '/Users/larour/ModelData/BedMap2/bedmap2_bin/'; +fid=fopen([path '/bedmap2_' string '.flt'],'r','l'); +data=fread(fid,[6667,6667],'float32'); +fclose(fid); + +% define grid +if strcmp(string,'thickness_uncertainty_5km'), + ncols =1361; + nrows =1361; + xll =-3401000; + yll =-3402000; + gridsize =5000; +else + ncols =6667; + nrows =6667; + xll =-3333000; + yll =-3333000; + gridsize =1000; +end +x_m=xll+(0:1:ncols-1)'*gridsize; +y_m=yll+(0:1:nrows-1)'*gridsize; + +%Change default to NaN +data(find(data==-9999))=NaN; + +%Interpolate +dataout = InterpFromGridToMesh(x_m,y_m,flipud(data'),double(X),double(Y),0); Index: /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/ismip6.m =================================================================== --- /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/ismip6.m (revision 27068) +++ /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/ismip6.m (revision 27069) @@ -17,5 +17,7 @@ thickness = {}; %placeholder for thicknesses deltathickness = {}; %placeholder for delta thicknesses + vdeltathickness = {}; %placeholder for delta thicknesses on vertices deltathicknessvaf = {}; %placeholder for delta thicknesses above floatation + vdeltathicknessvaf = {}; %placeholder for delta thicknesses above floatation on vertices deltathicknesshal = {}; %placeholder for delta thicknesses halosteric origins deltathicknessbar = {}; %placeholder for delta thicknesses halosteric origins @@ -80,4 +82,5 @@ fielddisplay(self,'thickness','thickness'); fielddisplay(self,'deltathickness','deltathickness'); + fielddisplay(self,'vdeltathickness','deltathickness at vertices'); fielddisplay(self,'icemask','icemask'); fielddisplay(self,'oceanmask','oceanmask'); Index: /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/oceanaverage.m =================================================================== --- /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/oceanaverage.m (revision 27069) +++ /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/oceanaverage.m (revision 27069) @@ -0,0 +1,18 @@ +function avg=oceanaverage(md,oceanlevelset,field) + + oceanel=isoceanin(oceanlevelset(md.mesh.elements)); + areas=GetAreasSphericalTria(md.mesh.elements,md.mesh.lat,md.mesh.long,planetradius('earth')); + pos=find(oceanel); + totalarea=sum(areas(pos)); + + if length(field)==md.mesh.numberofelements, + avg=sum(field(pos).*areas(pos))/totalarea; + elseif length(field)==md.mesh.numberofvertices, + avg=sum(field(md.mesh.elements(pos,:))*[1;1;1]/3.*areas(pos))/totalarea; + else error('field should be of size number of elements or numbrer of vertices'); + end + +end +function flag=isoceanin(ocean)% {{{ + flag=min(ocean,[],2)<0; +end% }}} Index: /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/oceanmip.m =================================================================== --- /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/oceanmip.m (revision 27068) +++ /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/oceanmip.m (revision 27069) @@ -14,7 +14,9 @@ zostoga = {}; % global-mean thermosteric sea level anomalies (in mm) (time) time = {}; %time in year + timem = {}; %time in months pbo = {}; %local ocean-bottom pressure changes. Also have zero mean over the oceans (lat,lon,time) lat = {}; long = {}; + mask = {}; scenario = {}; model = {}; @@ -44,7 +46,7 @@ %figure out time interval and remove historical: time=ncread(file,'time'); + pos=find(diff(time)<0); if isempty(pos), - %pos=(length(time)-12*100+1):length(time); pos=1:length(time); else @@ -55,4 +57,5 @@ pos=pos(pos2); time=time(pos); nt=length(time); + self.timem{end+1}=time; %reduce datasets: @@ -67,5 +70,5 @@ zosm(:,:,year)=mean(zos(:,:,(i-11):i),3); end - self.zos{end+1}=zosm; clear zos; + self.zos{end+1}=zosm; clear zos; clear zosm; %zostoga: @@ -92,4 +95,7 @@ self.long{end+1}=ncread(file,'lon'); + %mask: + self.mask{end+1}=ncread(file,'slm'); + %scenario and model: file=self.files{f}; @@ -102,5 +108,5 @@ end %}}} - function self = interpolate(self,md) % {{{ + function self = interpolate(self,md,slm) % {{{ %retrieve long and lat from mesh: @@ -110,4 +116,5 @@ for i=1:self.n, + %for i=1, disp(['interpolating model ' self.model{i} ' onto model mesh']); @@ -118,4 +125,5 @@ self.long{i}=self.long{i}(2:181,:,:); self.lat{i}=self.lat{i}(2:181,:,:); + self.mask{i}=self.mask{i}(2:181,:,:); end @@ -125,6 +133,12 @@ pos=find(long<0); long(pos)=long(pos)+360; lat=double(self.lat{i}); lat=lat(:); - [newl,uniqpos]=unique([lat,long],'rows','stable'); - long=long(uniqpos); lat=lat(uniqpos); + if slm==1, + mask=double(self.mask{i}); mask=mask(:); + maskpos=find(mask==1); + else + maskpos=(1:length(long))'; + end + [newl,uniqpos]=unique([lat(maskpos),long(maskpos)],'rows','stable'); + long=long(maskpos(uniqpos)); lat=lat(maskpos(uniqpos)); index=delaunay(long,lat); @@ -148,6 +162,6 @@ mesh_zos=zeros(md.mesh.numberofvertices,length(time)); - parfor j=1:length(time), - %for j=1:length(time), + %parfor j=1:length(time), + for j=1:length(time), if mod(j,10)==0, s=sprintf(' progress: %.2g ',j/length(time)*100); @@ -155,6 +169,6 @@ end - pbo=omip_pbo(:,:,j); pbo=pbo(:); pbo=pbo(uniqpos); - zos=omip_zos(:,:,j); zos=zos(:); zos=zos(uniqpos); + pbo=omip_pbo(:,:,j); pbo=pbo(:); pbo=pbo(maskpos(uniqpos)); + zos=omip_zos(:,:,j); zos=zos(:); zos=zos(maskpos(uniqpos)); pboj= InterpFromMeshToMesh2d(index,long,lat,pbo(newpos),meshlong,meshlat); Index: /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/solvebarystatic.m =================================================================== --- /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/solvebarystatic.m (revision 27069) +++ /issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/solvebarystatic.m (revision 27069) @@ -0,0 +1,99 @@ +function md=solvebarystatic(md) + + times=md.timestepping.start_time:md.timestepping.time_step:md.timestepping.final_time; + + nt=length(times); + SLR=zeros(nt,1); + VOL=zeros(nt,1); + totalarea=zeros(nt,1); + areas=GetAreasSphericalTria(md.mesh.elements,md.mesh.lat,md.mesh.long,planetradius('earth')); + + rho_ice=md.materials.rho_ice; + rho_water=md.materials.rho_water; + oceanarea=3.619e+14; + + responseice=zeros(md.mesh.numberofvertices,1); + responseocean=zeros(md.mesh.numberofvertices,1); + responsedh=zeros(md.mesh.numberofelements,nt-1); + responseeus=zeros(md.mesh.numberofelements,nt-1); + + for t=1:nt-1, + + time=times(t+1); + + %recover masks and thickness changes: + tdH=md.masstransport.spcthickness(end,:); + ind=findind(time,tdH); + dHtv=md.masstransport.spcthickness(1:end-1,ind)-md.geometry.thickness; + dHt=dHtv(md.mesh.elements)*[1;1;1]/3; + + tocean=md.mask.ocean_levelset(end,:); + ind=findind(time,tocean); + ocean=md.mask.ocean_levelset(1:end-1,ind); + + tice=md.mask.ice_levelset(end,:); + ind=findind(time,tice); + ice=md.mask.ice_levelset(1:end-1,ind); + + masks.isiceonly=isiceonly(ice(md.mesh.elements)); + masks.isoceanin=isoceanin(ocean(md.mesh.elements)); + masks.isfullyfloating=isfullyfloating(ocean(md.mesh.elements)); + masks.notfullygrounded=notfullygrounded(ocean(md.mesh.elements)); + + %Go through elements and sum up: + pos=find(masks.isiceonly & ~masks.isfullyfloating); + responsedh(:,t)=dHt; + responseice(:,t)=ice; + responseocean(:,t)=ocean; + responseeus(:,t)=ptve(md,pos); + SLR(t)=sum(-areas(pos).*dHt(pos)*rho_ice/(oceanarea*rho_water)); + + %total area: + totalarea(t)=sum(areas(pos)); + + %Volume too: + pos=find(masks.isiceonly); + VOL(t)=sum(-areas(pos).*dHt(pos)*rho_ice/(oceanarea*rho_water)); + %plotmodel(md,'data',masks.isiceonly,'view',[90 -90]); + %plotmodel(md,'data',areas.*dHt.*(masks.isiceonly & ~masks.isfullyfloating),'view',[90 -90]); + %plotmodel(md,'data',areas.*dHt.*(masks.isiceonly),'view',[90 -90]); + + end + + md.results.SLR=cumsum(SLR); + md.results.VOL=cumsum(VOL); + md.results.AREA=totalarea; + md.results.responseice=responseice; + md.results.responseocean=responseocean; + md.results.responsedh=responsedh; + md.results.responseeus=responseeus; +end +function ind=findind(t,ts) % {{{ + + nts=length(ts); + if t=ts(end), + ind=nts; + end; + for i=2:nts-1, + if (t>=ts(i) & t1, + tocean=md.mask.ocean_levelset(end,:); + ind=findind(time,tocean); + ocean=md.mask.ocean_levelset(1:end-1,ind); + else + ocean=md.mask.ocean_levelset; + end + + if size(md.mask.ice_levelset,2)>1, + tice=md.mask.ice_levelset(end,:); + ind=findind(time,tice); + ice=md.mask.ice_levelset(1:end-1,ind); + else + ice=md.mask.ice_levelset; + end + + masks.isiceonly=isiceonly(ice(md.mesh.elements)); + masks.isfullyfloating=isfullyfloating(ocean(md.mesh.elements)); + + %Go through elements and sum up: + %pos=find(~masks.isiceonly & ~masks.isfullyfloating); + %pos=1:md.mesh.numberofelements; + pos=find(~masks.isfullyfloating); + responsedh(:,t)=dHt; + responseice(:,t)=ice; + responseocean(:,t)=ocean; + responseeus(:,t)=ptve(md,pos); + SLR(t)=sum(-areas(pos).*dHt(pos)*rho_freshwater/(oceanarea*rho_water)); + + %total area: + totalarea(t)=sum(areas(pos)); + + end + md.results.SLR=cumsum(SLR); + md.results.AREA=totalarea; + md.results.responseice=responseice; + md.results.responseocean=responseocean; + md.results.responsedh=responsedh; + md.results.responseeus=responseeus; +end +function ind=findind(t,ts) % {{{ + + nts=length(ts); + if t=ts(end), + ind=nts; + end; + for i=2:nts-1, + if (t>=ts(i) & t