Index: /issm/trunk-jpl/src/m/model/averaging.m
===================================================================
--- /issm/trunk-jpl/src/m/model/averaging.m	(revision 12604)
+++ /issm/trunk-jpl/src/m/model/averaging.m	(revision 12605)
@@ -1,3 +1,3 @@
-function average=averaging(md,data,iterations)
+function average=averaging(md,data,iterations,varargin)
 %AVERAGING - smooths the input over the mesh
 %
@@ -8,33 +8,61 @@
 %   by taking the average of the element around a node weighted by the
 %   elements volume
+%   For 3d mesh, a last argument can be added to specify the layer to be averaged on.
 %
 %   Usage:
 %      smoothdata=averaging(md,data,iterations)
+%      smoothdata=averaging(md,data,iterations,layer)
 %
 %   Examples:
 %      velsmoothed=averaging(md,md.initialization.vel,4);
 %      pressure=averaging(md,md.initialization.pressure,0);
+%      temperature=averaging(md,md.initialization.temperature,1,1);
 
-if length(data)~=md.mesh.numberofelements & length(data)~=md.mesh.numberofvertices
+if ((nargin~=4) & (nargin~=3)),
+	error('averaging error message');
+end
+if (length(data)~=md.mesh.numberofelements & length(data)~=md.mesh.numberofvertices),
 	error('averaging error message: data not supported yet');
+end
+if md.mesh.dimension==3 & nargin==4,
+	if varargin{1}<=0 | varargin{1}>md.mesh.numberoflayers,
+		error('layer should be between 1 and md.mesh.numberoflayers');
+	end
+	layer=varargin{1};
+else
+	layer=0;
 end
 
 %initialization
-weights=zeros(md.mesh.numberofvertices,1);
-data=data(:);
+if layer==0,
+	weights=zeros(md.mesh.numberofvertices,1);
+	data=data(:);
+else 
+	weights=zeros(md.mesh.numberofvertices2d,1);
+	data=data((layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:);
+end
 
-%load some variables (it is much faster if the variab;es are loaded from md once for all)
-index=md.mesh.elements;
-numberofnodes=md.mesh.numberofvertices;
-numberofelements=md.mesh.numberofelements;
+%load some variables (it is much faster if the variabes are loaded from md once for all)
+if layer==0,
+	index=md.mesh.elements;
+	numberofnodes=md.mesh.numberofvertices;
+	numberofelements=md.mesh.numberofelements;
+else
+	index=md.mesh.elements2d;
+	numberofnodes=md.mesh.numberofvertices2d;
+	numberofelements=md.mesh.numberofelements2d;
+end
 
 %build some variables
 line=index(:);
-if md.mesh.dimension==3
+if md.mesh.dimension==3 & layer==0,
 	rep=6;
 	areas=GetAreas(index,md.mesh.x,md.mesh.y,md.mesh.z);
+elseif md.mesh.dimension==2,
+	rep=3;
+	areas=GetAreas(index,md.mesh.x,md.mesh.y);
 else
 	rep=3;
-	areas=GetAreas(index,md.mesh.x,md.mesh.y);
+	areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d);
 end
 summation=1/rep*ones(rep,1);