function [xc,yc,xg,yg,zA,DXG,DYG,DXC,DYC]=grid_sphere(dlon,dlat,varargin)
% Calculate grid data of fixed spherical grid
%
% [xc,yc,xg,yg,zA,DXG,DYG,DXC,DYC]=grid_sphere(dlon,dlat,lon0,lat0)

% Constants
Aearth=6370e3;

if nargin <=2
 lon0=0;
else
 lon0=varargin{1};
end
if nargin <= 3
 lat0=-90;
else
 lat0=varargin{2};
end
if prod(size(dlon)) == 1
 nx=round( 360/dlon );
 ddlon=dlon*ones([nx 1]);
else
 nx=prod(size(dlon));
 ddlon=reshape( dlon ,[nx 1]);
end
if prod(size(dlat)) == 1
 ny=round( (90-lat0)/dlat );
 ddlat=dlon*ones([1 ny]);
else
 ny=prod(size(dlat));
 ddlat=reshape( dlat ,[1 ny]);
end

ddlonc=[ddlon(1) (ddlon(1:end-1)'+ddlon(2:end)')/2 ddlon(end)]';
ddlatc=[ddlat(1) (ddlat(1:end-1)+ddlat(2:end))/2 ddlat(end)];

% Coordinates
xg=cumsum([lon0 ddlon']');
yg=cumsum([lat0 ddlat]);
xc=(xg(1:end-1)+xg(2:end))/2;
yc=(yg(1:end-1)+yg(2:end))/2;

if nargout <=4
 return
end

% Convert to radians
ddlon=ddlon*pi/180;
ddlat=ddlat*pi/180;
ddlonc=ddlonc*pi/180;
ddlatc=ddlatc*pi/180;

[XG,YG]=ndgrid(xg,yg);

clear XG
YU=(YG(:,1:end-1)+YG(:,2:end))/2;
YV=(YG(1:end-1,:)+YG(2:end,:))/2;
clear YG

% Surface area of grid-cell
zA=zeros(nx,ny);
for j=1:ny;
 zA(:,j)=Aearth^2*( ddlon.* ( sin( pi/180*YV(:,j+1) )-sin( pi/180*YV(:,j) ) ) );
end

if nargout <=5
 return
end

% Physical lengths only grid-cell edges
[DXG,DYC]=ndgrid(Aearth*ddlon,Aearth*ddlatc);
DXG=DXG.*cos(YV*pi/180);
[DXC,DYG]=ndgrid(Aearth*ddlonc,Aearth*ddlat);
DXC=DXC.*cos(YU*pi/180);

DXG( find(abs(DXG)<1.e-2) )=0;