function [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(nx,nratio,meth,ornt,plotfigs,writedata)

% Purser algorithm
% ngen=2; nrat2=1; nx=2*(1+nrat2*(2^ngen))-1;
% Number nodes along edge of cube
%nx=33;
% Ratio of working fine grid to target grid
%nratio=4;
% Method
%meth='conf';
%meth='q=1';
% Orientation
%ornt='c';
% Number of nodes along edge of fine cube (working grid)
nxf=nratio*(nx-1)+1;

%plotfigs=0;
%writedata=1;

% Generate gridded face using pure conformal mapping
%q=-1:2/(nxf-1):1;		% Coordinate -1 < q < +1 spans full face
 q=-1:2/(nxf-1):0;		% Coordinate -1 < q < 0 spans first quadrant
%q=rescale_coordinate(q,'q=0');
%q=rescale_coordinate(q,'q=78');
%q=rescale_coordinate(q,'q=1');
%q=rescale_coordinate(q,'q=i3');
%q=rescale_coordinate(q,'tan');
%q=rescale_coordinate(q,'conf');
 q=rescale_coordinate(q,meth);

[lx,ly]=ndgrid(q,q);		% Local fine grid curvilinear coordinate (lx,ly)

% Calculate simple finite volume info for fine grid (lx,ly);
[dx,dy,E] = calc_fvgrid(lx,ly);    clear dy;

% Use reflection symmetry of quadrants to fill face
dx((nxf+1)/2:nxf-1,:)=dx((nxf-1)/2:-1:1,:);
dx(:,(nxf+1)/2:nxf)=dx(:,(nxf+1)/2:-1:1);
E((nxf+1)/2:nxf-1,:)=E((nxf-1)/2:-1:1,:);
E(:,(nxf+1)/2:nxf-1)=E(:,(nxf-1)/2:-1:1);

[dxg,dxc,dxf,dxv]=reduce_dx(dx,nratio);
[Ec,Ez,Ev]=reduce_E(E,nratio);
clear dx E				% tidy up

% Remaining grid information
dyg=dxg';
dyc=dxc';
dyf=dxf';
dyu=dxv';
Eu=Ev';

% Calculate geographic coordinates
switch ornt
 case 'c'
  [LatG,LonG]=calc_geocoords_centerpole(lx,ly,1);
  for n=2:6; [LatG(:,:,n),LonG(:,:,n)]=calc_geocoords_centerpole(lx,ly,n); end;
  LatG=permutetiles(LatG,2);	% this is to be consistent with earlier grids
  LonG=permutetiles(LonG,2);	% this is to be consistent with earlier grids
 case 'v'
  [LatG,LonG]=calc_geocoords_cornerpole(lx,ly,1);
  for n=2:6; [LatG(:,:,n),LonG(:,:,n)]=calc_geocoords_cornerpole(lx,ly,n); end;
 otherwise
  ornt
  error('Unknown orientation');
end

if nratio~=1
% Calculate metric tensor
Q=q;Q(end:nxf)=-q(end:-1:1);
qx=Q(2+nratio/2:nratio:end)-Q(nratio/2:nratio:end-2);
[QX,QY]=ndgrid(qx,qx); clear qx Q
[Xf,Yf,Zf]=map_lonlat2xyz(LonG,LatG);
dXdx=( Xf(2+nratio/2:nratio:end  ,1+nratio/2:nratio:end,:) ...
      -Xf(  nratio/2:nratio:end-2,1+nratio/2:nratio:end,:) )./expand(QX);
dYdx=( Yf(2+nratio/2:nratio:end  ,1+nratio/2:nratio:end,:) ...
      -Yf(  nratio/2:nratio:end-2,1+nratio/2:nratio:end,:) )./expand(QX);
dZdx=( Zf(2+nratio/2:nratio:end  ,1+nratio/2:nratio:end,:) ...
      -Zf(  nratio/2:nratio:end-2,1+nratio/2:nratio:end,:) )./expand(QX);
dXdy=( Xf(1+nratio/2:nratio:end,2+nratio/2:nratio:end  ,:) ...
      -Xf(1+nratio/2:nratio:end,  nratio/2:nratio:end-2,:) )./expand(QY);
dYdy=( Yf(1+nratio/2:nratio:end,2+nratio/2:nratio:end  ,:) ...
      -Yf(1+nratio/2:nratio:end,  nratio/2:nratio:end-2,:) )./expand(QY);
dZdy=( Zf(1+nratio/2:nratio:end,2+nratio/2:nratio:end  ,:) ...
      -Zf(1+nratio/2:nratio:end,  nratio/2:nratio:end-2,:) )./expand(QY);
Q11=dXdx.*dXdx+dYdx.*dYdx+dZdx.*dZdx;
Q22=dXdy.*dXdy+dYdy.*dYdy+dZdy.*dZdy;
Q12=dXdx.*dXdy+dYdx.*dYdy+dZdx.*dZdy;
clear Xf Yf Zf QX QY Q
else
Q11=0;
Q12=0;
Q22=0;
end

% Sub-sample to obtain model coordinates
latG=LatG(1:nratio:end,1:nratio:end,:);
lonG=LonG(1:nratio:end,1:nratio:end,:);
if nratio==1
 latC=( latG(1:end-1,1:end-1,:) ...
       +latG(2:end  ,1:end-1,:) ...
       +latG(1:end-1,2:end  ,:) ...
       +latG(2:end  ,2:end  ,:) )/4;
 lonC=( lonG(1:end-1,1:end-1,:) ...
       +lonG(2:end  ,1:end-1,:) ...
       +lonG(1:end-1,2:end  ,:) ...
       +lonG(2:end  ,2:end  ,:) )/4;
else
 latC=LatG(1+nratio/2:nratio:end,1+nratio/2:nratio:end,:);
 lonC=LonG(1+nratio/2:nratio:end,1+nratio/2:nratio:end,:);
end
clear LatG LonG				% tidy up

if nratio~=1
% Flow rotation tensor
Xlon=-sin(lonC); Ylon= cos(lonC); Zlon= 0*lonC;
TUu= (dXdx.*Xlon+dYdx.*Ylon+dZdx.*Zlon);
TVu=-(dXdy.*Xlon+dYdy.*Ylon+dZdy.*Zlon);
Xlat=-sin(latC).*cos(lonC); Ylat=-sin(latC).*sin(lonC); Zlat= cos(latC);
TUv=-(dXdx.*Xlat+dYdx.*Ylat+dZdx.*Zlat);
TVv= (dXdy.*Xlat+dYdy.*Ylat+dZdy.*Zlat);
det=sqrt(TUu.*TVv-TUv.*TVu);
%det=sqrt(TUu.*TUu+TUv.*TUv);
TUu=TUu./det;
TUv=TUv./det;
%det=sqrt(TVv.*TVv+TVu.*TVu);
TVu=TVu./det;
TVv=TVv./det;
clear dXdx dXdy dYdx dYdy dZdx dZdy
end

% 3D coordinates for plotting
[XG,YG,ZG]=map_lonlat2xyz(lonG,latG);

% Symmetry?
stats( dxg-dxg(end:-1:1,:) )
stats( dxg(:,1:end)-dxg(:,end:-1:1) )
stats( dyg-dyg(:,end:-1:1) )
stats( dyg(1:end,:)-dyg(end:-1:1,:) )
stats( dxg-dyg' )
stats( dxf-dxf(end:-1:1,:) )
stats( dxf-dxf(:,end:-1:1) )
stats( dyf-dyf(end:-1:1,:) )
stats( dxf-dyf' )
stats( Ec-Ec(end:-1:1,:) )
stats( Ec-Ec(:,end:-1:1) )
stats( Ec-Ec' )
stats( Eu-Ev' )
stats( Eu-Eu(:,end:-1:1) )
stats( Ev-Ev(end:-1:1,:) )

if writedata~=0
dir=sprintf('cs-%s-%s-%i-%i/',meth,ornt,nx-1,nxf-1);
[status,msg]=mkdir(dir); if status==0; disp(msg); end;
clear status msg

fid=fopen([dir 'grid.info'],'w');
fprintf(fid,'nx=%i\n',nx-1);
fprintf(fid,'nratio=%i\n',nratio);
fprintf(fid,'nxf=%i\n',(nx-1)*nratio);
fprintf(fid,'mapping="%s"\n',meth);
fprintf(fid,'orientation="%s"\n',ornt);
fclose(fid);

Rsphere=6370e3;
write_tiles([dir 'LONG'], 180/pi*lonG )
write_tiles([dir 'LATG'], 180/pi*latG )
write_tiles([dir 'LONC'], 180/pi*lonC )
write_tiles([dir 'LATC'], 180/pi*latC )
%write_tiles([dir 'DXG'], Rsphere*expand(dxg) )
%write_tiles([dir 'DYG'], Rsphere*expand(dyg) )
%write_tiles([dir 'DXF'], Rsphere*expand(dxf) )
%write_tiles([dir 'DYF'], Rsphere*expand(dyf) )
%write_tiles([dir 'DXC'], Rsphere*expand(dxc) )
%write_tiles([dir 'DYC'], Rsphere*expand(dyc) )
%write_tiles([dir 'DXV'], Rsphere*expand(dxv) )
%write_tiles([dir 'DYU'], Rsphere*expand(dyu) )
%write_tiles([dir 'RA'], Rsphere^2*expand(Ec) )
%write_tiles([dir 'RAW'], Rsphere^2*expand(Eu) )
%write_tiles([dir 'RAS'], Rsphere^2*expand(Ev) )
%write_tiles([dir 'RAZ'], Rsphere^2*expand(Ez) )
write_tile([dir 'DXG'], Rsphere*(dxg) )
write_tile([dir 'DYG'], Rsphere*(dyg) )
write_tile([dir 'DXF'], Rsphere*(dxf) )
write_tile([dir 'DYF'], Rsphere*(dyf) )
write_tile([dir 'DXC'], Rsphere*(dxc) )
write_tile([dir 'DYC'], Rsphere*(dyc) )
write_tile([dir 'DXV'], Rsphere*(dxv) )
write_tile([dir 'DYU'], Rsphere*(dyu) )
write_tile([dir 'RA'], Rsphere^2*(Ec) )
write_tile([dir 'RAW'], Rsphere^2*(Eu) )
write_tile([dir 'RAS'], Rsphere^2*(Ev) )
write_tile([dir 'RAZ'], Rsphere^2*(Ez) )

% Save stuff for pre-processing/post-processing
save([dir 'CUBE_HGRID.mat'],'latG','lonG','latC','lonC','dxg','dyg','dxc','dyc','Ec');
save([dir 'CUBE_3DGRID.mat'],'XG','YG','ZG');
save([dir 'TUV.mat'],'TUu','TUv','TVu','TVv');
end

if plotfigs~=0
% Some basic pictures
figure(1);
subplot(221);imagesc(dxv');axis xy;title('\Delta x_v');colorbar
subplot(222);imagesc(dxg');axis xy;title('\Delta x_g');colorbar
subplot(223);imagesc(dxc');axis xy;title('\Delta x_c');colorbar
subplot(224);imagesc(dxf');axis xy;title('\Delta x_f');colorbar

figure(2);
subplot(221);imagesc( ((Ec-dxf.*dxf')./Ec)' );axis xy;title('|E_c-\Delta x_f\Delta y_f|');colorbar
subplot(222);imagesc(Ec');axis xy;title('E_c');colorbar
subplot(223);imagesc(Ez');axis xy;title('E_z');colorbar
subplot(224);imagesc(Ev');axis xy;title('E_v');colorbar

% Plot rotation tensor
figure(3)
subplot(221)
plotcube(XG,YG,ZG,TUu);colorbar;%[az el]=view;view([az+180 el])
title('TUu')
subplot(222)
plotcube(XG,YG,ZG,TUv);colorbar;%[az el]=view;view([az+180 el])
title('TUv')
subplot(223)
plotcube(XG,YG,ZG,TVu);colorbar;%[az el]=view;view([az+180 el])
title('TVu')
subplot(224)
plotcube(XG,YG,ZG,TVv);colorbar;%[az el]=view;view([az+180 el])
title('TVv')
end