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	cnh	1.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;