format compact
path('/u/u0/dimitri/matlab/tools',path);
path('/u/u0/dimitri/matlab/whoi',path);
path('/u/u0/dimitri/matlab/mitgcm',path);
path('/u/adcroft/matlab/bin',path);

nx=32;

%for typ={'conf','q=0','q=1','q=1/2','q=7/8','q=i3','tan','tan2','new'}

for typ={'conf'}

[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+1,16,typ{1},'c',0,1);
theta=zeros(nx/2,nx/2);
for i=1:nx/2, for j=1:nx/2
  lat=[latG(i,j) latG(i+1,j)]*180/pi;
  lon=[lonG(i,j) lonG(i+1,j)]*180/pi;
  a=dist(lat,lon,'sphere');
  lat=[latG(i,j) latG(i,j+1)]*180/pi;
  lon=[lonG(i,j) lonG(i,j+1)]*180/pi;
  b=dist(lat,lon,'sphere');
  lat=[latG(i+1,j) latG(i,j+1)]*180/pi;
  lon=[lonG(i+1,j) lonG(i,j+1)]*180/pi;
  c=dist(lat,lon,'sphere');
  theta(i,j)=acos((a^2+b^2-c^2)/2/a/b);
end, end
disp(typ{1})
disp(minmax(dxg*6371))
disp(minmax(dxg(1:nx/2,1:nx/2)./dyg(1:nx/2,1:nx/2)))
tmp=theta*180/pi; tmp(1,1)=nan;
disp(minmax(tmp))
mypcolor(theta*180/pi); colorbar
title(typ{1})
tmp=typ{1}; tmp(strfind(tmp,'/'))='_';
eval(['print -dpsc FIG' tmp])

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all, close all, nx=32; genx;
pn='/u/u0/dimitri/MITgcm/verification/hs94.cs-32x32x5/input/';
PN='/u/u0/dimitri/cubed_sphere/cs-conf-c-32-128/';

long=readbin([pn 'LONG.bin'],[nx 6 nx],1,'real*8');
latg=readbin([pn 'LATG.bin'],[nx 6 nx],1,'real*8');
LONG=zeros(nx,6,nx); LATG=zeros(nx,6,nx);
for i=1:6
 tmp=readbin([PN 'LONG.' myint2str(i,3) '.bin'],[nx+1 nx+1],1,'real*8');
 LONG(:,i,:)=tmp(1:nx,1:nx);
 tmp=readbin([PN 'LATG.' myint2str(i,3) '.bin'],[nx+1 nx+1],1,'real*8');
 LATG(:,i,:)=tmp(1:nx,1:nx);
end
for i=1:6
 subplot(3,2,i)
 plot(squeeze(long(:,i,:)),squeeze(latg(:,i,:)),'bo', ...
      squeeze(LONG(:,i,:)),squeeze(LATG(:,i,:)),'k.')
 title(i)
end

clf
subplot(211), fld=[pn 'LONG.bin']; plotfld, title('LONG')
subplot(212), fld=[pn 'LATG.bin']; plotfld, title('LATG')

lonc=readbin([pn 'LONC.bin'],[nx 6 nx],1,'real*8');
latc=readbin([pn 'LATC.bin'],[nx 6 nx],1,'real*8');
LONC=zeros(nx,6,nx); LATC=zeros(nx,6,nx);
for i=1:6
 LONC(:,i,:)=readbin([PN 'LONC.' myint2str(i,3) '.bin'],[nx nx],1,'real*8');
 LATC(:,i,:)=readbin([PN 'LATC.' myint2str(i,3) '.bin'],[nx nx],1,'real*8');
end
for i=1:6
 subplot(3,2,i)
 plot(squeeze(lonc(:,i,:)),squeeze(latc(:,i,:)),'bo', ...
      squeeze(LONC(:,i,:)),squeeze(LATC(:,i,:)),'k.')
 title(i)
end

for nm={'DXC','DYC','DXF','DYF','DXG','DYG','DXV','DYU','RA','RAS','RAW','RAZ'}
 fld=tile_fld(PN,nm{1},nx);
 writebin([PN nm{1} '.bin'],fld,1,'real*8');
 subplot(311), fld=[pn nm{1} '.bin']; plotfld, title([pn nm{1}]), tmp1=tmp;
 subplot(312), fld=[PN nm{1} '.bin']; plotfld, title([PN nm{1}]), tmp2=tmp;
 subplot(313), mypcolor(tmp2'-tmp1'); colorbar, title('difference')
 disp(minmax((tmp1-tmp2)./tmp2)), pause
end