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

%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 degrees'), end;
if max(lon)<=2*pi; error('Error- longitude input should be in degress'), end; 


if length(find(lon==0))>1 & (find(lon==0)) > 1  %longitude passes through zero
    k=find(lon==0);
    lon=[lon(1:k-1);lon(k:N)+360];
end
    
    
out=zeros(NN,N);

a=6.37*10^6;

%convert to radians
lon=lon*pi/180;
lat=lat*pi/180;
[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)= (-1*vc(3,m)+4*vc(2,m)-3*vc(1,m))/(2*(lat(2)-lat(1)));
    vy(NN,m)= (1*vc(NN-2,m)-4*vc(NN-1,m)+3*vc(NN,m))/(2*(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));

for m=1:1:NN
    ux(m,1)= (-1*U(m,3)+4*U(m,2)-3*U(m,1))/(2*(a*cos(lat(m))*lon(2)-a*cos(lat(m))*lon(1)));
    ux(m,N)= (1*U(m,N-2)-4*U(m,N-1)+3*U(m,N))/(2*(a*cos(lat(m))*lon(N)-a*cos(lat(m))*lon(N-1)));
end
    
for m=2:1:N-1
    for n=1:1:NN
        ux(n,m)=(U(n,m+1)-U(n, m-1))/(a*cos(lat(n))*(lon(m+1)-lon(m-1)));
    end
end

out=ux+vy;
   