Index: /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/confidenceintervals.m
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/confidenceintervals.m	(revision 25588)
+++ /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/confidenceintervals.m	(revision 25588)
@@ -0,0 +1,36 @@
+function level=confidenceintervals(x,h,threshold)
+%CONFIDENCEINTERVALS: build confidence interval from histogram.
+   
+	nods=size(h,1); nbins=size(h,2);
+
+	%make sure h is one column smaller than x
+	if size(x,2) ~= (size(h,2)+1),
+		error('x should be one column larger than the histogram');
+	end
+
+	%First build cdf: 
+	c=[zeros(nods,1) cumsum(h,2)];
+
+	%In percentage: 
+	threshold=threshold/100;
+
+	flags=zeros(nods,nbins+1);
+	pos=find(c<threshold);
+	flags(pos)=1;
+	diffc=diff(flags,1,2);
+
+	[ind,j]=find(diffc==-1);
+	indices=zeros(nods,1); 
+	indices(ind)=j;
+
+	%threshold was passed between j and j+1 for each node. try 
+	%to pin down the fraction better: 
+	idx = sub2ind(size(c), [1:size(c,1)]', indices);
+	idx2 = sub2ind(size(c), [1:size(c,1)]', indices+1);
+	val1=c(idx);
+	val2=c(idx2);
+	x1=x(idx);
+	x2=x(idx2);
+
+	fraction=(threshold-val1)./(val2-val1);
+	level=x1+fraction.*(x2-x1);
Index: /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/glacier_inventory.m	(revision 25587)
+++ /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/glacier_inventory.m	(revision 25588)
@@ -149,4 +149,5 @@
 			errornotfound=getfieldvalue(options,'errornotfound',1);
 			plotr=getfieldvalue(options,'plot',0);
+			plotr=1;
 
 			%The way we run this is through the O2 zones defined in boxes. We go through 
@@ -203,8 +204,12 @@
 					switch i,
 						case 4, radius=40;
+						case 8, radius=60;
 						case 12, radius=25;
+						case 19, radius=60;
+						case 32, radius=60;
 						case 33, radius=5;
 						case 41, radius=75;
 						case 42, radius=45;
+						case 61, radius=66;
 						case 68, radius=10;
 						case 82, radius=30;
Index: /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/oceanmip.m
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/oceanmip.m	(revision 25587)
+++ /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/oceanmip.m	(revision 25588)
@@ -46,5 +46,6 @@
 					pos=find(diff(time)<0); 
 					if isempty(pos),
-						pos=(length(time)-12*100+1):length(time);
+						%pos=(length(time)-12*100+1):length(time);
+						pos=1:length(time);
 					else
 						pos=[(pos+1):length(time)]; 
@@ -75,4 +76,6 @@
 						zostogam(year)=mean(zostoga(i-11:i));
 					end
+					%control against 2006: 
+					zostogam=zostogam-zostogam(1);
 					self.zostoga{end+1}=zostogam; clear zostoga;
 
@@ -109,6 +112,16 @@
 				disp(['interpolating model ' self.model{i} ' onto model mesh']);
 
+				%If we have 182 cells in long, trim by 1 on each side:
+				if size(self.pbo{i},1)==182,
+					self.pbo{i}=self.pbo{i}(2:181,:,:);
+					self.zos{i}=self.zos{i}(2:181,:,:);
+					self.long{i}=self.long{i}(2:181,:,:);
+					self.lat{i}=self.lat{i}(2:181,:,:);
+				end
+
+
 				%triangulation: 
-				long=double(self.long{i});  long=long(:);
+				long=double(self.long{i});  long=long(:); 
+				pos=find(long<0); long(pos)=long(pos)+360;
 				lat=double(self.lat{i}); lat=lat(:);
 				[newl,uniqpos]=unique([lat,long],'rows','stable');
@@ -122,4 +135,7 @@
 				ind=find(areas<1e-8); index(ind,:)=[];
 
+				%fix if we have orphans
+				[index long lat dummy newpos]=FixMesh(index,long,lat,1:length(long));
+
 				time=self.time{i};
 				%retrieve fields:
@@ -133,4 +149,5 @@
 
 				parfor j=1:length(time),
+				%for j=1:length(time),
 					if mod(j,10)==0, 
 						s=sprintf('   progress: %.2g ',j/length(time)*100);
@@ -138,9 +155,14 @@
 					end
 
-					pbo=omip_pbo(:,:,j); pbo=pbo(:); pbo=pbo(uniqpos);
+					pbo=omip_pbo(:,:,j); pbo=pbo(:); pbo=pbo(uniqpos); 
 					zos=omip_zos(:,:,j); zos=zos(:); zos=zos(uniqpos);
 
-					mesh_pbo(:,j)=InterpFromMeshToMesh2d(index,long,lat,pbo,meshlong,meshlat);
-					mesh_zos(:,j)=InterpFromMeshToMesh2d(index,long,lat,zos,meshlong,meshlat);
+					pboj= InterpFromMeshToMesh2d(index,long,lat,pbo(newpos),meshlong,meshlat);
+					pos=find(abs(pboj)>1e4); pboj(pos)=0;
+					mesh_pbo(:,j)=pboj;
+
+					zosj=InterpFromMeshToMesh2d(index,long,lat,zos(newpos),meshlong,meshlat);
+					pos=find(abs(zosj)>1e5); zosj(pos)=0;
+					mesh_zos(:,j)=zosj;
 				end
 				self.mesh_pbo{end+1}=mesh_pbo;
@@ -200,6 +222,5 @@
 						pbo=self.mesh_pbo{i}; pbo=pbo/1000; %in meters
 					end
-					time=self.time{i};
-					array=[pbo;time];
+					array=pbo;
 					break;
 				end