code obtained from faulks.lcs.mit.edu:/scratch/adcroft/cubed_sphere

path('/home/dimitri/matlab/adcroft/bin',path);
tic
[dxg,dyg,dxf,dyf,dxc,dyc,dxv,dyu,Ec,Eu,Ev,Ez,latC,lonC,latG,lonG,...
          Q11,Q22,Q12, TUu,TUv,TVu,TVv ]=gengrid_fn(576,4,'conf','c',0,0);
toc

timings                           nireas  faulks
gengrid_fn(576,2,'conf','c',0,0)  47.9    60.0
gengrid_fn(576,4,'conf','c',0,0)  163.9   

current 80S-80N, 1/4-deg isotropic grid:
1440 * 1088 = 1566720 grid points
27.7987 km at Equator, 4.8272 km at 80N

tripolar grid:
2304 grid boxes around Equator
2*pi*6371/2304 = 17.3742 km grid at Equator
switch to bipolar is at 60N where grid is 17.3742*cos(60) =  8.6871 km
at North Pole grid size is 17.3742/pi = 5.5304 km
smallest distance in Canadian Archipelago is approximately 17.3742/6 = 2.8957 km
worst aspect ratio is approximately 3
number of grid points in sperical polar region is 2304*960 = 2211840
number of grid points in bipolar grid is 1152*1152 = 1327104
total number of grid points is 3538944
advantages, grid is perfectly orthogonal
there are no singularities
no modifications to sea-ice thermodynamics

by comparison
the cubed sphere with 2304 grid boxes around Equator
using Alistair's eq 16 has:
576*576*6 = 1990656 grid points (56% of tripolar grid, 127% of 1/4-deg grid)
largest distance is 22.9683 (132% of tripolar grid)
smallest distance is 3.2633 (112% of tripolar grid, 68% of 1/4-deg grid)
worst aspect ratio is 1.6162 (53% of tripolar grid)
but orthogonality is not as good as that of tripolar grid
and sea-ice dynamics need some work
also for the tripolar grid, removing land point will
result in much more dramatic reduction in number of grid
points because the three singularities are over land.



table
first row is min, max, and mean of dxg and dyg
second row is min, max, and mean of dxg./dyg
third row is min, max, and mean of angles, excluding cube corners

conf, nratio=4
    2.1460   18.0293   15.9182
    0.9172    1.0903    1.0000
   89.9978   98.5344   90.0562

conf, nratio=2
    2.1460   18.0293   15.9182
    0.9167    1.0909    1.0000
   89.9978   98.5344   90.0562

q=0, nratio=4
    2.2544   17.5336   15.9487
    0.9172    1.0903    1.0006
   89.9977   98.5344   90.0554

q=0, nratio=2
    2.2544   17.4069   15.9468
    0.9167    1.0909    1.0006
   89.9977   98.5344   90.0554

q=1, nratio=4 (has NaNs)
   13.6379   67.1398   16.1959
    0.2031    4.9226    1.0552
   89.9972   98.5343   90.0497

q=1, nratio=2
   13.6310   67.1393   16.1925
    0.2031    4.9238    1.0551
   89.9972   98.5343   90.0497

q=1/2, nratio=4
    2.8990   21.0171   16.0604
    0.6992    1.4302    1.0137
   89.9974   98.5344   90.0527

q=1/2, nratio=2
    2.8958   21.0001   16.0613
    0.6999    1.4287    1.0136
   89.9974   98.5344   90.0527

q=7/8, nratio=4 (has NaNs)
   13.6379   67.1398   16.1959
    0.2031    4.9226    1.0552
   89.9972   98.5343   90.0497

q=7/8, nratio=2
   13.6310   67.1393   16.1925
    0.2031    4.9238    1.0551
   89.9972   98.5343   90.0497

q=i3, nratio=4
   13.1144   65.1937   16.1921
    0.2092    4.7807    1.0549
   89.9972   98.5343   90.0497

q=i3, nratio=2
   12.4259   62.6094   16.1915
    0.2179    4.5899    1.0546
   89.9973   98.5343   90.0498

tan, nratio=4
    3.2633   22.9683   16.0912
    0.6188    1.6162    1.0203
   89.9974   98.5344   90.0519

tan, nratio=2
    3.2633   22.9683   16.0912
    0.6188    1.6162    1.0203
   89.9974   98.5344   90.0519

tan2, nratio=4
    2.2237   17.7883   15.9327
    0.9174    1.0900    1.0001
   89.9978   98.5344   90.0558

tan2, nratio=2
    2.2237   17.7883   15.9327
    0.9169    1.0906    1.0001
   89.9978   98.5344   90.0558

new, nratio=4
    2.5095   18.8618   15.9191
    0.8889    1.1250    1.0000
   89.9978   98.5344   90.0561

new, nratio=2
    2.8863   20.9484   15.9201
    0.8000    1.2500    1.0002
   89.9978   98.5344   90.0561

 




gengrids					% Does it all
gengrids_fn(nx,nratio,projection,orientation,plotfigs,writedata)
[dxg,dyg,dxf,dyf,dxc,dyc,dxv,dyu,Ec,Eu,Ev,Ez,latC,lonC,latG,lonG,...
          Q11,Q22,Q12, TUu,TUv,TVu,TVv ]=gengrid_fn(24,2,'conf','c',1,0);

[X,Y,Z,x,y]=read_scubdat(fn,n);			% Read SCUB0002.DAT file
[X,Y,Z]=map_xy2xyz(x,y);			% Map (x,y) -> (X,Y,Z)
[X,Y,Z]=gmap_xy2xyz(x,y);			% Map (x,y) -> (X,Y,Z)

[lon,lat]=map_xyz2lonlat(X,Y,Z);		% Map (X,Y,Z) -> (lon,lat)
[X,Y,Z]=map_lonlat2xyz(lon,lat);		% Map (lon,lat) -> (X,Y,Z)
[lat,lon]=calc_geocoords(lx,ly,tile)		% Map (x,y) -> (lat,lon)

[X,Y,Z]=rotate_about_xaxis(lX,lY,lZ,angle);	% Rotate about X axis
[X,Y,Z]=rotate_about_yaxis(lX,lY,lZ,angle);	% Rotate about Y axis
[X,Y,Z]=rotate_about_zaxis(lX,lY,lZ,angle);	% Rotate about Z axis

V=coord2vector(XG,YG,ZG);			% Convert [X,Y,Z] to vector V

angle=angle_between_coords(X1,Y1,Z1,X2,Y2,Z2);	% Angles between coords
angle=angle_between_vectors(V1,V2);		% Angles between vectors V1,V2
excess=excess_of_quad(V1,V2,V3,V4);		% Excess of quadrilateral

q=rescale_coordinate(q,'q=i3');			% Re-scale single coordinate
[dx,dy,E]=calc_fvgrid(lx,ly);			% Calculate F.V. grid info
dxg=conf_dx(q);					% Calculate conformal dxg
[dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio);		% Sum fine-grid dx -> dxg,...
[Ec,Ez,Ev]=reduce_E(dx,nratio);			% Sum fine-grid E -> Ec,...

a=bari(b);					% Two point average in I dir.
a=barj(b);					% Two point average in J dir.

write_tiles('LATC',lat)				% Write tiled files
displaytiles(lat)				% Display tiles in comp. layout
plotcube(X,Y,Z,C)				% Display tiles as 3-D sphere
merccube(lonG,latG,C)				% Display tiles Mercator proj.