function [urot,vrot]=rotate(lat,lon,um,vm,u,v)
%function [urot,vrot]=rotate(lat,lon,um,vm,u,v)
%
% --- Rotate vectors represented on a general grid to meridional and
% --- zonal componets.
%
% --- Created 2004 by mats.bentsen@nersc.no
%

nx=size(lat,1);
ny=size(lat,2);          

rad=1.7453292519943295E-02;

urot=zeros(nx,ny);
vrot=zeros(nx,ny); 

for j=1:ny
  jp1=mod(j,ny)+1;
  for i=1:nx
    ii=max([2,min([nx-1,i])]);
    dlat=lat(ii+1,j)-lat(ii-1,j);
    dlon=mod(lon(ii+1,j)+360.,360.) ...
      -mod(lon(ii-1,j)+360.,360.);
    if (abs(dlon+360.)<abs(dlon)) 
      dlon=dlon+360.;
    end
    if (abs(dlon-360.)<abs(dlon))
      dlon=dlon-360.;
    end
    phi=atan2(dlat,cos(lat(ii,j)*rad)*dlon); 
    cosphi=cos(phi);
    sinphi=sin(phi);
    ip1=min([nx,i+1]);
    up=(u(i,j)*um(i,j)+u(ip1,j)*um(ip1,j)) ...
       /max([1.,um(i,j)+um(ip1,j)]);
    vp=(v(i,j)*vm(i,j)+v(i,jp1)*vm(i,jp1)) ...
       /max([1.,vm(i,j)+vm(i,jp1)]);
    urot(i,j)=up*cosphi-vp*sinphi;
    vrot(i,j)=up*sinphi+vp*cosphi;
  end
end

