function [int]=zon_int(array,X,Y)
%zonal integration of a 2D geophysical field using trapezoidal rule
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%array : the array to work on, of dimension ny*nx
%X : longitude in degrees (vector (nx))
%Y : latitude in degrees (vector (ny))
%lat : the latitude at which the integration is performed, in degrees.
%returns int(Y).
%RMK : works for Y scalar as well.
%JEG, July 2002

conv=pi/180.;
R=6371000.;
nx=size(X,1);
ny=size(Y,1);

if (ny > 1)  
	for jj=1:ny
		integrd=cos(conv*Y(jj))*array(jj,:);
		int(jj)=R*conv*trapz(X,integrd);
	end
else 	integrd=cos(conv*Y)*array(:);
		int=R*conv*trapz(X,integrd);
end