function [out]=div(U,V,lon,lat)

%calculates the divergence of a vector field in spherical geometry 
c=size(U);
N=c(1,2);
NN=c(1,1);


if nargin<4; error('Error - Not enough input parameters specified'), end;
if nargin>4; error('Error - to0 many input parameters specified'), end;
if size(U)~=size(V); error('Error- vector component matrices must be the same size'), end;
if length(lat)>NN;error('Error - vector of latitude values has too large of dimensions'),end;
if length(lat)<NN;error('Error - vector of latitude values has too small of dimensions'),end;
if length(lon)>N;error('Error - vector of longitude values has too large of dimensions'),end;
if length(lon)<NN;error('Error - vector of longitude values has too small of dimensions'),end;
if max(lat)>pi; error('Error- latitude input should be in radians'), end;
if max(lon)>2*pi; error('Error- longitude input should be in radians'), end; 

out=zeros(NN,N);

a=6.37*10^6;

[LON LAT]=meshgrid(lon,lat);
vc=V.*cos(LAT);  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%IMPORTANT CHNAGE FROM RECTANGULAR%%%%%%%%%%%

vy=zeros(NN,N);

for m=1:1:N
    vy(1,m)= (vc(2,m)-vc(1,m))/(lat(2)-lat(1));
    vy(NN,m)= (vc(NN,m)-vc(NN-1,m))/(lat(NN)-lat(NN-1));
end
    
for m=2:1:NN-1
    for n=1:1:N
        vy(m,n)=(vc(m+1,n)-vc(m-1,n))/(lat(m+1)-lat(m-1));
    end
end

vy=vy./(a*cos(LAT));

ux=divx(U, lon, lat);

out=ux+vy;
   