Index: ../trunk-jpl/src/m/contrib/larour/glaciermip.m
===================================================================
--- ../trunk-jpl/src/m/contrib/larour/glaciermip.m	(nonexistent)
+++ ../trunk-jpl/src/m/contrib/larour/glaciermip.m	(revision 24848)
@@ -0,0 +1,152 @@
+%GLACIERMIP class definition
+%
+%   Usage:
+%      gmip = glaciermip('ncfile',glaciermipncfile,'version',1)
+%
+%
+
+classdef glaciermip < handle
+	properties (SetAccess=public) %Model fields
+		ncfile  = '';
+		version = NaN;
+
+		%version 2:
+		area =[];
+		mass =[];
+		region =[];
+		time =[];
+		glaciermodel =[];
+		climatemodel =[];
+		scenario =[];
+
+		%from version 1:
+		run=[];
+		forcingmodel=[];
+		realization=[];
+		volume=[];
+	end
+	methods
+		function self = glaciermip(varargin) % {{{
+
+			if nargin==0, 
+				self=setdefaultparameters(self);
+			else 
+				self=setdefaultparameters(self);
+
+				options=pairoptions(varargin{:});
+
+				self.ncfile=getfieldvalue(options,'ncfile');
+				self.version=getfieldvalue(options,'version');
+
+
+				%read variables:
+				if self.version==1,
+					self.region=ncread(self.ncfile,'region');
+					self.time=ncread(self.ncfile,'time');
+					self.run=ncread(self.ncfile,'run');
+					self.glaciermodel=ncread(self.ncfile,'glaciermodel');
+					self.forcingmodel=ncread(self.ncfile,'forcingmodel');
+					self.scenario=ncread(self.ncfile,'scenario');
+					self.realization=ncread(self.ncfile,'realization');
+					self.area=ncread(self.ncfile,'area');
+					self.volume=ncread(self.ncfile,'volume');
+				elseif self.version==2,
+					self.area=ncread(self.ncfile,'Area');
+					self.mass=ncread(self.ncfile,'Mass');
+					self.region=ncread(self.ncfile,'Region');
+					self.time=ncread(self.ncfile,'Time');
+					self.glaciermodel=ncread(self.ncfile,'Glacier_Model');
+					self.climatemodel=ncread(self.ncfile,'Climate_Model');
+					self.scenario=ncread(self.ncfile,'Scenario');
+					%mass(region,time,climatemodel,glaciermodel,scenario)
+				else 
+					error(sprintf('glaciermipfile constructor error message: version %i for MIP not supported!'),self.version);
+				end
+			end
+		end
+		%}}}
+		function inv = setdefaultparameters(inv) % {{{
+		end
+		%}}}
+		function t=gettime(self,varargin) % {{{
+			
+			options=pairoptions(varargin{:});
+			t=self.time;
+
+		end % }}}
+		function masses=getmass(self,varargin) % {{{
+			
+			options=pairoptions(varargin{:});
+
+			rg=getfieldvalue(options,'region',1:length(self.region));
+			cm=getfieldvalue(options,'climatemodel',1:length(self.climatemodel));
+			gm=getfieldvalue(options,'glaciermodel',1:length(self.glaciermodel));
+			sc=getfieldvalue(options,'scenario',1:length(self.scenario));
+			zerostonan=getfieldvalue(options,'zerostonan',0);
+
+			if self.version==1,
+				error(sprintf('getmass not supported yet for Glacier MIP version %i',self.version));
+			end
+
+			%serialize: 
+			masses=[];
+			for i=rg,
+				for j=cm,
+					for k=gm,
+						for l=sc,
+							masses=[masses; squeeze(self.mass(i,:,j,k,l))];
+						end
+					end
+				end
+			end
+			if zerostonan,
+				masses(find(masses==0))=NaN;
+			end
+
+		end % }}}
+		function massrates=getmassrates(self,varargin) % {{{
+		
+			options=pairoptions(varargin{:});
+
+			%get mass first: 
+			masses=self.getmass(varargin{:});
+
+			%compute mass rates: 
+			dt=diff(self.time);
+			dm=diff(masses,1,2);
+
+			massrates=dm;
+			for i=1:size(massrates,1),
+				massrates(i,:)= massrates(i,:)./dt';
+			end
+
+
+		end % }}}
+		function disp(self) % {{{
+			disp(sprintf('   Glacier MIP (version %i):',self.version)); 
+
+			if self.version==1,
+				fielddisplay(self,'ncfile','netcdf file for GlacierMIP results');
+				fielddisplay(self,'time','time scale in yr');
+				fielddisplay(self,'region','region');
+				fielddisplay(self,'run','run');
+				fielddisplay(self,'glaciermodel','glacier model');
+				fielddisplay(self,'forcingmodel','forcing model');
+				fielddisplay(self,'scenario','scenario');
+				fielddisplay(self,'realization','realization');
+				fielddisplay(self,'area','area');
+				fielddisplay(self,'volume','volume');
+			end 
+			if self.version==2,
+				fielddisplay(self,'ncfile','netcdf file for GlacierMIP results');
+				fielddisplay(self,'time','time');
+				fielddisplay(self,'region','region');
+				fielddisplay(self,'glaciermodel','glaciermodel');
+				fielddisplay(self,'climatemodel','climatemodel');
+				fielddisplay(self,'scenario','scenario');
+				fielddisplay(self,'area','area');
+				fielddisplay(self,'mass','mass');
+			end
+		end % }}}
+	end
+end
Index: ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24847)
+++ ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24848)
@@ -92,6 +92,22 @@
 			end
 		end
 		%}}}
+		function totalarea=area(self,varargin) % {{{
+			region=-1;
+			totalarea=0;
+			if nargin==2,
+				region=varargin{1};
+			end
+			if region==-1,
+				%figure out the areas of everybody: 
+				for i=1:self.nregions(),
+					totalarea=totalarea+sum(self.regions(i).Area);
+				end
+			else
+				totalarea=totalarea+sum(self.regions(region).Area);
+			end
+
+		end % }}}
 		function [mpartition,npartition]=partition(self,varargin) % {{{
 			mpartition=zeros(self.nregions(),1);
 			npartition=zeros(self.nregions(),1);
@@ -239,7 +255,7 @@
 		%}}}
 		function mesh_connectivity(self,mesh) % {{{
 			
-			%initialize glacier_connectivity: 
+			%initialize glacier and element connectivity: 
 			glacier_connectivity=zeros(self.nglaciers,1);
 
 			%assume maximum width for connectivity table and initialize 
@@ -325,14 +341,11 @@
 			%For each element, this will tell us which glaciers belong to this element. 
 			for j=1:length(glacier_connectivity),
 				el=glacier_connectivity(j);
-				if el==0,
-					continue;
-				end
 				count=element_connectivity(el,ny);
 				if count>ny,
 					error('need to enlarge connectivity table');
 				end
-				element_connectivity(el,count+1)=partition_counter+j;
+				element_connectivity(el,count+1)=j;
 				element_connectivity(el,ny)=count+1;
 			end