function meanfield = mit_zonalmean(field,hfac,dx);
%function meanfield = mit_zonalmean(field,hfac,dx);

% $Header: /home/ubuntu/mnt/e9_copy/MITgcm/verification/tutorial_global_oce_latlon/diags_matlab/mit_zonalmean.m,v 1.3 2006/08/12 20:25:13 jmc Exp $
% $Name: checkpoint65d $

  [nx ny nz] = size(field);

  area = hfac.*repmat(dx,[1 1 nz]);
  meanfield = squeeze(nanmean(field.*area)./nanmean(area));
  
  knil = find(meanfield == 0);
  meanfield(knil) = NaN;
  
  return