function [Tm] = zonalmean(GRID,T,varargin);
% Tm = zonalmean(GRID,T,varargin);

T( find(isnan(T)) )=0;

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

Tm=squeeze(sum(dTm.*repmat(mapmsk,[1 1 size(T,3)]),1));
N=squeeze(sum( GRID.mskc.*repmat(mapmsk,[1 1 size(T,3)]),1));
Tm=Tm./sq(N);