Index: /issm/trunk-jpl/src/m/contrib/larour/oceanmip.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/oceanmip.m	(revision 25001)
+++ /issm/trunk-jpl/src/m/contrib/larour/oceanmip.m	(revision 25001)
@@ -0,0 +1,197 @@
+%OCEANMIP class definition
+%
+%   Usage:
+%      omip = oceanmip('root',rootdir,'files',listfiles);
+%
+
+classdef oceanmip < handle
+	properties (SetAccess=public) %Model fields
+		
+		root = ''; %where are the files for CMIP5
+		n   = 0;   %number of files
+		files   = {};   %netcdf files
+		zos = {}; % local sea-level anomalies with zero mean over the oceans (in mm) (lat,lon,time): 
+		zostoga = {}; % global-mean thermosteric sea level anomalies (in mm) (time)
+		time = {}; %time in year
+		pbo = {}; %local ocean-bottom pressure changes. Also have zero mean over the oceans (lat,lon,time)
+		lat = {};
+		long = {};
+		scenario = {}; 
+		model = {}; 
+		mesh_zos = {}; %interpolated (onto model mesh) zos field
+		mesh_pbo = {}; %interpolated (onto model mesh) pbo field
+
+	end
+	methods
+		function self = oceanmip(varargin) % {{{
+
+			if nargin==0, 
+				self=setdefaultparameters(self);
+			else 
+				self=setdefaultparameters(self);
+
+				options=pairoptions(varargin{:});
+
+				self.root=getfieldvalue(options,'root');
+				self.files=getfieldvalue(options,'files');
+				self.n=length(self.files);
+
+				%read files: 
+				for f=1:self.n,
+					file=[self.root '/' self.files{f}];
+					disp(['reading file ' file]);
+					
+					%figure out time interval: 
+					time=ncread(file,'time');
+					pos=find(time>2000); time=time(pos); nt=length(time);
+					
+					%reduce datasets: 
+					time=floor(time(12:12:nt)); 
+					self.time{end+1}=time;
+					
+					%zos: 
+					zos=ncread(file,'zos'); nx=size(zos,1); ny=size(zos,2); zos=zos(:,:,pos);
+					zosm=zeros(nx,ny,nt/12);
+					for i=12:12:nt,
+						year=i/12;
+						zosm(:,:,year)=mean(zos(:,:,i-11:i),3);
+					end
+					self.zos{end+1}=zosm; clear zos;
+
+					%zostoga:
+					zostoga=ncread(file,'zostoga');
+					zostogam=zeros(nt/12);
+					for i=12:12:nt,
+						year=i/12;
+						zostogam(year)=mean(zostoga(i-11:i));
+					end
+					self.zostoga{end+1}=zostogam; clear zostoga;
+
+					%pbo: 
+					pbo=ncread(file,'pbo'); nx=size(pbo,1); ny=size(pbo,2); pbo=pbo(:,:,pos);
+					pbom=zeros(nx,ny,nt/12);
+					for i=12:12:nt,
+						year=i/12;
+						pbom(:,:,year)=mean(pbo(:,:,i-11:i),3);
+					end
+					self.pbo{end+1}=pbom; clear pbo;
+
+					self.lat{end+1}=ncread(file,'lat');
+					self.long{end+1}=ncread(file,'lon');
+
+					%scenario and model: 
+					file=self.files{f};
+					ind=findstr(file,'rcp'); 
+					self.scenario{end+1}=str2num(file(ind+3:ind+4));
+					self.model{end+1}=file(1:ind-2);
+				end
+
+			end
+		end
+		%}}}
+		function self = interpolate(self,md) % {{{
+
+			%retrieve long and lat from mesh: 
+			meshlong=md.mesh.long;
+			pos=find(meshlong<0); meshlong(pos)=meshlong(pos)+360;
+			meshlat=md.mesh.lat;
+
+			for i=1:self.n,
+				disp(['interpolating model ' self.model{i} ' onto model mesh']);
+
+				%triangulation: 
+				long=double(self.long{i});  long=long(:);
+				lat=double(self.lat{i}); lat=lat(:);
+				[newl,uniqpos]=unique([lat,long],'rows','stable');
+				long=long(uniqpos); lat=lat(uniqpos);
+				index=delaunay(long,lat);
+
+				%some checks: 
+				areas=GetAreas(index,long,lat); indneg=find(areas<0);  
+				index(indneg,:)=[index(indneg,1),index(indneg,3),index(indneg,2)];
+				areas=GetAreas(index,long,lat); 
+				ind=find(areas<1e-8); index(ind,:)=[];
+
+				time=self.time{i};
+				%retrieve fields:
+				omip_pbo=self.pbo{i};
+				omip_zos=self.zos{i};
+				omip_zostoga=self.zostoga{i};
+
+				%interpolate:
+				mesh_pbo=zeros(md.mesh.numberofvertices,length(time));
+				mesh_zos=zeros(md.mesh.numberofvertices,length(time));
+
+				parfor j=1:length(time),
+					if mod(j,10)==0, 
+						s=sprintf('   progress: %.2g ',j/length(time)*100);
+						fprintf(s); pause(1e-3); fprintf(repmat('\b',1,numel(s))); 
+					end
+
+					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);
+				end
+				self.mesh_pbo{end+1}=mesh_pbo;
+				self.mesh_zos{end+1}=mesh_zos;
+
+				%clear for memory purposes: 
+				self.pbo{i}=[];
+				self.zos{i}=[];
+			end
+		end
+
+		%}}}
+		function self = setdefaultparameters(self) % {{{
+		end
+		%}}}
+		function self = rawclean(self) % {{{
+			for i=1:self.n,
+				self.zos{i}=[];
+				self.pbo{i}=[];
+			end
+		end
+		%}}}
+		function [pbo,time]= bottompressure(self,model,gridded) % {{{
+			for i=1:self.n,
+				if strcmpi(model,self.model{i}),
+					if gridded,
+						pbo=self.pbo{i}; pbo=pbo/1000; %in meters
+					else
+						pbo=self.mesh_pbo{i}; pbo=pbo/1000; %in meters
+					end
+					time=self.time{i};
+					break;
+				end
+			end
+		end
+		%}}}
+		function [lat,long]= latlong(self,model) % {{{
+			for i=1:self.n,
+				if strcmpi(model,self.model{i}),
+					lat=self.lat{i};
+					long=self.long{i};
+					break;
+				end
+			end
+		end
+		%}}}
+		function disp(self) % {{{
+			disp('   CMIP5 Ocean MIP:');
+
+				fielddisplay(self,'n','number of files');
+				fielddisplay(self,'root','where are the files for CMIP5');
+				fielddisplay(self,'files','CMIP5 ocean files');
+				fielddisplay(self,'zos','local sea-level anomalies with zero mean over the oceans (in mm) (lat,lon,time)');
+				fielddisplay(self,'zostoga','global-mean thermosteric sea level anomalies (in mm) (time)');
+				fielddisplay(self,'time','time in years');
+				fielddisplay(self,'pbo','local ocean-bottom pressure changes. Also have zero mean over the oceans (lat,lon,time)');
+				fielddisplay(self,'lat','latitudes');
+				fielddisplay(self,'long','longitudes');
+				fielddisplay(self,'scenario','scenarios');
+				fielddisplay(self,'model','model names');
+		end % }}}
+	end 
+end