function [Tm,Tsd] = meanprofile(GRID,T1,varargin);

% Remove NaN's
T1( find(isnan(T1)) )=0;

if nargin==3
 dTm=(mean(T1,4)-mean(varargin{1},4)).*GRID.mskc;
 mapmsk=GRID.mskc(:,:,1);
elseif nargin==4
 dTm=(mean(T1,4)-mean(varargin{1},4)).*GRID.mskc;
 mapmsk=varargin{2};
else
 dTm=mean(T1,4).*GRID.mskc;
 mapmsk=GRID.mskc(:,:,1);
end

nxy=prod(size(GRID.rac));
nr=prod(size(GRID.mskc,3));
ra=reshape((GRID.rac(:).*mapmsk(:))*ones(1,nr),size(GRID.mskc)).*GRID.mskc;
dTm=reshape(dTm,[nxy nr]);
ra=reshape(ra,[nxy nr]);
Tm=squeeze( sum(dTm.*ra) ./ sum(ra) );
Tsd=(dTm'-Tm'*ones(1,nxy))';
Tsd=sqrt( squeeze( sum((Tsd.^2).*ra) ./ sum(ra) ) );

function [F] = impliedflux(GRID,Q);
F=mean(Q,3).*GRID.mskc(:,:,1).*GRID.rac;
F=-cumsum([0 squeeze( sum(F,1) )]);