Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24843)
+++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24843)
@@ -0,0 +1,122 @@
+%GLACIER_INVENTORY class definition
+%
+%   Usage:
+%      inv = glacier_inventory(varargin)
+%
+%      where varargin is a variable list of options: 
+%
+%   Example: 
+%      rgi = glacier_inventory('shapefileroot',shapefileroot,'region_names',region_names);
+%
+%
+
+classdef glacier_inventory < handle
+	properties (SetAccess=public) %Model fields
+		region_names     = {};
+		shapefileroot    = '';
+		regions          = struct();
+	end
+	methods
+		
+		function self = glacier_inventory(varargin) % {{{
+
+			if nargin==0, 
+				self=setdefaultparameters(self);
+			else 
+				self=setdefaultparameters(self);
+
+				options=pairoptions(varargin{:}); 
+
+				self.shapefileroot=getfieldvalue(options,'shapefileroot');
+				self.region_names=getfieldvalue(options,'region_names');
+
+				%read the shape files and create the regions: 
+				counter=0;
+				self.regions=struct();
+				for i=1:self.nregions(),
+					self.regions(i).name=self.region_names{i};
+					disp(['reading region: '  self.regions(i).name]);
+					self.regions(i).id=i;
+					contours=shpread([self.shapefileroot '/' self.regions(i).name '.shp']);
+					
+					%we can't keep all the info: build arrays of centroids:
+					O1Region=zeros(length(contours),1);
+					O2Region=zeros(length(contours),1);
+					Area=zeros(length(contours),1);
+					CenLon=zeros(length(contours),1);
+					CenLat=zeros(length(contours),1);
+					for j=1:length(contours),
+						O1Region(j)=str2num(contours(j).O1Region);
+						O2Region(j)=str2num(contours(j).O2Region);
+						Area(j)=contours(j).Area;
+						CenLon(j)=contours(j).CenLon;
+						CenLat(j)=contours(j).CenLat;
+					end
+
+					self.regions(i).Area=Area;
+					self.regions(i).O1Region=O1Region;
+					self.regions(i).O2Region=O2Region;
+					self.regions(i).CenLat=CenLat;
+					self.regions(i).CenLon=CenLon;
+					self.regions(i).lids=[1:length(contours)]';
+					self.regions(i).gids=self.regions(i).lids+counter;
+					counter=counter+length(contours);
+				end
+			end
+		end
+		%}}}
+		function inv = setdefaultparameters(inv) % {{{
+		end
+		%}}}
+		function n = nregions(self) % {{{
+			n=length(self.region_names);
+		end
+		%}}}
+		function disp(self) % {{{
+			disp(sprintf('   Glacier inventory:')); 
+
+			disp(['   number of regions: ' num2str(self.nregions())]);
+			for i=1:self.nregions(),
+				disp(sprintf('      region %i: ''%s'' %i glaciers ',i,self.regions(i).name,length(self.regions(i).Area)));
+			end
+
+		end % }}}
+		function checkconsistency(slm,solutiontype) % {{{
+
+			%is the coupler turned on? 
+			for i=1:length(slm.icecaps),
+				if slm.icecaps{i}.transient.iscoupler==0,
+					warning(sprintf('sealevelmodel checkconsistenty error:  icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));
+				end
+			end
+				
+			if slm.earth.transient.iscoupler==0,
+				warning('sealevelmodel checkconsistenty error:  earth model should have the transient coupler option turned on!');
+			end
+
+			%check that the transition vectors have the right size: 
+			for i=1:length(slm.icecaps),
+				if slm.icecaps{i}.mesh.numberofvertices ~= length(slm.earth.slr.transitions{i}),
+					error(['sealevelmodel checkconsistenty issue with size of transition vector for ice cap: ' num2str(i) ' name: ' slm.icecaps{i}.miscellaneous.name]);
+				end
+			end
+			
+			%check that run_frequency is the same everywhere: 
+			for i=1:length(slm.icecaps),
+				if slm.icecaps{i}.slr.geodetic_run_frequency~=slm.earth.slr.geodetic_run_frequency,
+					error(sprintf('sealevelmodel checkconsistenty error:  icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name));
+				end
+			end
+
+			%make sure steric_rate is the same everywhere: 
+			for i=1:length(slm.icecaps),
+				md= slm.icecaps{i}; 
+				if ~isempty(find(md.slr.steric_rate - slm.earth.slr.steric_rate(slm.earth.slr.transitions{i}))),
+					error(sprintf('steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
+				end
+			end
+
+		end
+		%}}}
+	end
+end