% generate rotated  u/v vectors and compare to model output
% and to those used by global_ocean.cs32x15

clear all, close all

% load model grid
xc=readbin('XC.data',[32,6,32]);xc=permute(xc,[1 3 2]);
xg=readbin('XG.data',[32,6,32]);xg=permute(xg,[1 3 2]);
yc=readbin('YC.data',[32,6,32]);yc=permute(yc,[1 3 2]);
yg=readbin('YG.data',[32,6,32]);yg=permute(yg,[1 3 2]);
%tmp=yg; cx=[min(tmp(:)) max(tmp(:))]; plot_cube

% load input velocity on Gaussian NCEP grid
u=readbin('uwind_192_94_12.bin',[192 94]);
v=readbin('vwind_192_94_12.bin',[192 94]);
lon=0:1.875:358.125;
lat=[-88.5420, -86.6532, -84.7532, -82.8508, -80.9474, -79.0435, ...
     -77.1393, -75.2351, -73.3307, -71.4262, -69.5217, -67.6171, ...
     -65.7125, -63.8079, -61.9033, -59.9986, -58.0940, -56.1893, ...
     -54.2846, -52.3799, -50.4752, -48.5705, -46.6658, -44.7611, ...
     -42.8564, -40.9517, -39.0470, -37.1422, -35.2375, -33.3328, ...
     -31.4281, -29.5234, -27.6186, -25.7139, -23.8092, -21.9044, ...
     -19.9997, -18.0950, -16.1902, -14.2855, -12.3808, -10.4760, ...
     - 8.5713, - 6.6666, - 4.7618, - 2.8571, - 0.9524, ...
      0.9524, 2.8571, 4.7618, 6.6666, 8.5713, 10.4760, 12.3808, ...
      14.2855, 16.1902, 18.0950, 19.9997, 21.9044, 23.8092, ...
      25.7139, 27.6186, 29.5234, 31.4281, 33.3328, 35.2375, ...
      37.1422, 39.0470, 40.9517, 42.8564, 44.7611, 46.6658, ...
      48.5705, 50.4752, 52.3799, 54.2846, 56.1893, 58.0940, ...
      59.9986, 61.9033, 63.8079, 65.7125, 67.6171, 69.5217, ...
      71.4262, 73.3307, 75.2351, 77.1393, 79.0435, 80.9474, ...
      82.8508, 84.7532, 86.6532, 88.5420];

% interpolate u and v on xc, yc
u1=[u(192,:); u; u(1,:)];
u1=[u1(:,1) u1 u1(:,94)];
v1=[v(192,:); v; v(1,:)];
v1=[v1(:,1) v1 v1(:,94)];
lon1=-1.875:1.875:360;
lat1=[-90 lat 90];
u2=interp2(lon1,lat1,u1',xc(:),yc(:));
v2=interp2(lon1,lat1,v1',xc(:),yc(:));
u2=reshape(u2,size(xc)); v2=reshape(v2,size(xc));

% compute adjacent grid coordinates for dx
nx=1:31; ny=1:31;
x1=xg(nx,ny,:);
x2=xg(nx+1,ny,:);
x3=xg(nx,ny+1,:);
x4=xg(nx+1,ny+1,:);
ix=find(x2-x1>180); x2(ix)=x2(ix)-360;
ix=find(x1-x2>180); x2(ix)=x2(ix)+360;
ix=find(x3-x1>180); x3(ix)=x3(ix)-360;
ix=find(x1-x3>180); x2(ix)=x3(ix)+360;
ix=find(x4-x1>180); x4(ix)=x4(ix)-360;
ix=find(x1-x4>180); x4(ix)=x4(ix)+360;

% compute adjacent grid coordinates for dy
y1=yg(nx,ny,:);
y2=yg(nx+1,ny,:);
y3=yg(nx,ny+1,:);
y4=yg(nx+1,ny+1,:);

% compute cube-sphere u
dx=0.5*(x2+x4-x1-x3).*cos(yc(nx,ny,:)*pi/180); minmax(dx)
dy=0.5*(y2+y4-y1-y3); minmax(dy)
denom=1./sqrt(dx.*dx+dy.*dy);
u3=(u2(nx,ny,:).*dx+v2(nx,ny,:).*dy).*denom;

% compute cube-sphere v
dx=0.5*(x3+x4-x1-x2).*cos(yc(nx,ny,:)*pi/180); minmax(dx)
dy=0.5*(y3+y4-y1-y2); minmax(dy)
denom=1./sqrt(dx.*dx+dy.*dy);
v3=(u2(nx,ny,:).*dx+v2(nx,ny,:).*dy).*denom;

cx=[-10 10];
figure(1), tmp=u3; plot_cube, title('u from matlab')
figure(2), tmp=v3; plot_cube, title('v from matlab')

% load global_ocean.cs32x15 input files
un=permute(readbin('../inp_thsice/ncep_uwind_cs.bin',[32,6,32],1,'real*8'),[1 3 2]);
figure(3), tmp=un; cx=[min(tmp(:)) max(tmp(:))]; plot_cube
title('u from ocean.cs32x15')
print -djpeg u_cs32
vn=permute(readbin('../inp_thsice/ncep_vwind_cs.bin',[32,6,32],1,'real*8'),[1 3 2]);
figure(4), tmp=vn; cx=[min(tmp(:)) max(tmp(:))]; plot_cube
title('v from ocean.cs32x15')
print -djpeg v_cs32

% load model output files
um=permute(readbin('UWIND.0000000020.data',[32,6,32]),[1 3 2]);
figure(5), tmp=um; plot_cube, title('u from model output')
print -djpeg u_exf
vm=permute(readbin('VWIND.0000000020.data',[32,6,32]),[1 3 2]);
figure(6), tmp=vm; plot_cube, title('v from model output')
print -djpeg v_exf