function [mn]=idma_area_mean(fld,msk);
% IDMA_AREA_MEAN computes area-weighted average of a
%    field (fld; 2D or 3D) over a region defined by 
%    a 0/1 mask (msk). Both fld and msk are of 
%    the gcmfaces class. If fld has a third dimension
%    (e.g., for depth or time) then mn is a vector
%
% Example: box_l=[-180 -140]; box_L=[-10 10]-30;
%          msk=(mygrid.XC>=box_l(1)&mygrid.XC<box_l(2)&...
%               mygrid.YC>=box_L(1)&mygrid.YC<box_L(2));
%
%          fileName='release1/nctiles_climatology/THETA/THETA'; 
%          fldName='THETA';
%          fld=read_nctiles(fileName,fldName,1);
%          fld=fld.*mygrid.mskC;
%
%          [mn]=idma_area_mean(fld,msk);
%          figureL; plot(mn,mygrid.RC);

gcmfaces_global;

n=size(fld{1},3);
mn=NaN*zeros(1,n);

racmsk=msk.*mygrid.RAC;
racmsk=repmat(racmsk,[1 1 n]);
racmsk(isnan(fld))=0;

for k=1:n;
    mn(k)=nansum(racmsk(:,:,k).*fld(:,:,k))/nansum(racmsk(:,:,k));
end;