| 1 | 
enderton | 
1.1 | 
function [AngleCS,AngleSN] = cubeCalcAngle(YG,RAC,dxG,dyG); | 
| 2 | 
  | 
  | 
% [AngleCS,AngleSN] = cubeCalcAngle(YG,RAC,dxG,dyG); | 
| 3 | 
  | 
  | 
% | 
| 4 | 
  | 
  | 
% Determine rotation angle, alpha (returned as cos(alpha) and sin(alpha)), | 
| 5 | 
  | 
  | 
% to convert cube sphere C-grid vectors to east-west and north-south | 
| 6 | 
  | 
  | 
% components. | 
| 7 | 
  | 
  | 
% | 
| 8 | 
  | 
  | 
% e.g. | 
| 9 | 
  | 
  | 
% | 
| 10 | 
  | 
  | 
% >> XG=rdmds('XG'); | 
| 11 | 
  | 
  | 
% >> YG=rdmds('YG'); | 
| 12 | 
  | 
  | 
% >> RAC=rdmds('RAC'); | 
| 13 | 
  | 
  | 
% >> dxG=rdmds('dxG'); | 
| 14 | 
  | 
  | 
% >> dyG=rdmds('dyG'); | 
| 15 | 
  | 
  | 
% >> [AngleCS,AngleSN] = cubeCalcAngle(YG,RAC,dxG,dyG); | 
| 16 | 
  | 
  | 
% | 
| 17 | 
  | 
  | 
% >> u=rdmds('uVeltave.0000513360'); | 
| 18 | 
  | 
  | 
% >> v=rdmds('vVeltave.0000513360'); | 
| 19 | 
  | 
  | 
% >> [uE,vN] = rotate_csCg_EN(u,v,AngleCS,AngleSN,XG,YG); | 
| 20 | 
  | 
  | 
 | 
| 21 | 
  | 
  | 
nc=size(RAC,2); | 
| 22 | 
  | 
  | 
deg2rad=pi/180; | 
| 23 | 
  | 
  | 
r_earth=sqrt(sum(sum(RAC))./(4.*pi)); | 
| 24 | 
  | 
  | 
 | 
| 25 | 
  | 
  | 
[yZ6]=split_Z_cub(YG); | 
| 26 | 
  | 
  | 
 | 
| 27 | 
  | 
  | 
% Build purely zonal flow from a stream function, psi = r_earth * latitude. | 
| 28 | 
  | 
  | 
psi=-r_earth.*yZ6.*deg2rad; | 
| 29 | 
  | 
  | 
uZ=psi([1:nc],[1:nc],:)-psi([1:nc],[2:nc+1],:); | 
| 30 | 
  | 
  | 
vZ=psi([2:nc+1],[1:nc],:)-psi([1:nc],[1:nc],:); | 
| 31 | 
  | 
  | 
uZ=reshape(permute(uZ,[1 3 2]),[6.*nc nc])./dyG; | 
| 32 | 
  | 
  | 
vZ=reshape(permute(vZ,[1 3 2]),[6.*nc nc])./dxG; | 
| 33 | 
  | 
  | 
 | 
| 34 | 
  | 
  | 
% Put constructed zonal wind at cell center. | 
| 35 | 
  | 
  | 
[uu,vv]=split_UV_cub(uZ,vZ); | 
| 36 | 
  | 
  | 
u=(uu(1:nc,:,:)+uu(2:nc+1,:,:))/2; | 
| 37 | 
  | 
  | 
v=(vv(:,1:nc,:)+vv(:,2:nc+1,:))/2; | 
| 38 | 
  | 
  | 
uZc=reshape(permute(u,[1 3 2]),[6.*nc nc]); | 
| 39 | 
  | 
  | 
vZc=reshape(permute(v,[1 3 2]),[6.*nc nc]); | 
| 40 | 
  | 
  | 
 | 
| 41 | 
  | 
  | 
% Calculate angles. | 
| 42 | 
  | 
  | 
norm=sqrt(uZc.*uZc+vZc.*vZc); | 
| 43 | 
  | 
  | 
AngleCS =  uZc./norm; | 
| 44 | 
  | 
  | 
AngleSN = -vZc./norm; |