matlab-grid-generator/bin/grid_sphere.m

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;
1	function [xc,yc,xg,yg,zA,DXG,DYG,DXC,DYC]=grid_sphere(dlon,dlat,varargin)
2	% Calculate grid data of fixed spherical grid
3	%
4	% [xc,yc,xg,yg,zA,DXG,DYG,DXC,DYC]=grid_sphere(dlon,dlat,lon0,lat0)
5
6	% Constants
7	Aearth=6370e3;
8
9	if nargin <=2
10	lon0=0;
11	else
12	lon0=varargin{1};
13	end
14	if nargin <= 3
15	lat0=-90;
16	else
17	lat0=varargin{2};
18	end
19	if prod(size(dlon)) == 1
20	nx=round( 360/dlon );
21	ddlon=dlon*ones([nx 1]);
22	else
23	nx=prod(size(dlon));
24	ddlon=reshape( dlon ,[nx 1]);
25	end
26	if prod(size(dlat)) == 1
27	ny=round( (90-lat0)/dlat );
28	ddlat=dlon*ones([1 ny]);
29	else
30	ny=prod(size(dlat));
31	ddlat=reshape( dlat ,[1 ny]);
32	end
33
34	ddlonc=[ddlon(1) (ddlon(1:end-1)'+ddlon(2:end)')/2 ddlon(end)]';
35	ddlatc=[ddlat(1) (ddlat(1:end-1)+ddlat(2:end))/2 ddlat(end)];
36
37	% Coordinates
38	xg=cumsum([lon0 ddlon']');
39	yg=cumsum([lat0 ddlat]);
40	xc=(xg(1:end-1)+xg(2:end))/2;
41	yc=(yg(1:end-1)+yg(2:end))/2;
42
43	if nargout <=4
44	return
45	end
46
47	% Convert to radians
48	ddlon=ddlon*pi/180;
49	ddlat=ddlat*pi/180;
50	ddlonc=ddlonc*pi/180;
51	ddlatc=ddlatc*pi/180;
52
53	[XG,YG]=ndgrid(xg,yg);
54
55	clear XG
56	YU=(YG(:,1:end-1)+YG(:,2:end))/2;
57	YV=(YG(1:end-1,:)+YG(2:end,:))/2;
58	clear YG
59
60	% Surface area of grid-cell
61	zA=zeros(nx,ny);
62	for j=1:ny;
63	zA(:,j)=Aearth^2( ddlon. ( sin( pi/180YV(:,j+1) )-sin( pi/180YV(:,j) ) ) );
64	end
65
66	if nargout <=5
67	return
68	end
69
70	% Physical lengths only grid-cell edges
71	[DXG,DYC]=ndgrid(Aearthddlon,Aearthddlatc);
72	DXG=DXG.cos(YVpi/180);
73	[DXC,DYG]=ndgrid(Aearthddlonc,Aearthddlat);
74	DXC=DXC.cos(YUpi/180);
75
76	DXG( find(abs(DXG)<1.e-2) )=0;