Changeset 27069

Timestamp:

06/16/22 17:14:58 (3 years ago)

Author:

Eric.Larour

Message:

CHG: diverse.

Location:

issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour

Files:

• interpBedmap2.m (added)
• ismip6.m (modified) (2 diffs)
• oceanaverage.m (added)
• oceanmip.m (modified) (11 diffs)
• solvebarystatic.m (added)
• solvebarystatic2.m (added)
• solvebarystatichydro.m (added)

issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/ismip6.m

@@ -17 +17 @@
                 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 +82 @@
                         fielddisplay(self,'thickness','thickness');
                         fielddisplay(self,'deltathickness','deltathickness');
+                        fielddisplay(self,'vdeltathickness','deltathickness at vertices');
                         fielddisplay(self,'icemask','icemask');
                         fielddisplay(self,'oceanmask','oceanmask');

issm/branches/trunk-larour-SLPS2022/src/m/contrib/larour/oceanmip.m

@@ -14 +14 @@
                 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 +46 @@
                                         %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 +57 @@
                                         pos=pos(pos2);
                                         time=time(pos); nt=length(time);
+                                        self.timem{end+1}=time;
 
                                         %reduce datasets: 
@@ -67 +70 @@
                                                 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 +95 @@
                                         self.long{end+1}=ncread(file,'lon');
 
+                                        %mask: 
+                                        self.mask{end+1}=ncread(file,'slm'); 
+
                                         %scenario and model: 
                                         file=self.files{f};
@@ -102 +108 @@
                 end
                 %}}}
-                function self = interpolate(self,md) % {{{
+                function self = interpolate(self,md,slm) % {{{
 
                         %retrieve long and lat from mesh: 
@@ -110 +116 @@
 
                         for i=1:self.n,
+                        %for i=1,
                                 disp(['interpolating model ' self.model{i} ' onto model mesh']);
 
@@ -118 +125 @@
                                         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 +133 @@
                                 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 +162 @@
                                 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 +169 @@
                                         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);