Diagnostics/DiagUtility/cubeCalcAngle.m

function [AngleCS,AngleSN] = cubeCalcAngle(YG,RAC,dxG,dyG);
% [AngleCS,AngleSN] = cubeCalcAngle(YG,RAC,dxG,dyG);
%
% Determine rotation angle, alpha (returned as cos(alpha) and sin(alpha)),
% to convert cube sphere C-grid vectors to east-west and north-south
% components.
%
% e.g.
%
% >> XG=rdmds('XG');
% >> YG=rdmds('YG');
% >> RAC=rdmds('RAC');
% >> dxG=rdmds('dxG');
% >> dyG=rdmds('dyG');
% >> [AngleCS,AngleSN] = cubeCalcAngle(YG,RAC,dxG,dyG);
%
% >> u=rdmds('uVeltave.0000513360');
% >> v=rdmds('vVeltave.0000513360');
% >> [uE,vN] = rotate_csCg_EN(u,v,AngleCS,AngleSN,XG,YG);

nc=size(RAC,2);
deg2rad=pi/180;
r_earth=sqrt(sum(sum(RAC))./(4.*pi));

[yZ6]=split_Z_cub(YG);

% Build purely zonal flow from a stream function, psi = r_earth * latitude.
psi=-r_earth.*yZ6.*deg2rad;
uZ=psi([1:nc],[1:nc],:)-psi([1:nc],[2:nc+1],:);
vZ=psi([2:nc+1],[1:nc],:)-psi([1:nc],[1:nc],:);
uZ=reshape(permute(uZ,[1 3 2]),[6.*nc nc])./dyG;
vZ=reshape(permute(vZ,[1 3 2]),[6.*nc nc])./dxG;

% Put constructed zonal wind at cell center.
[uu,vv]=split_UV_cub(uZ,vZ);
u=(uu(1:nc,:,:)+uu(2:nc+1,:,:))/2;
v=(vv(:,1:nc,:)+vv(:,2:nc+1,:))/2;
uZc=reshape(permute(u,[1 3 2]),[6.*nc nc]);
vZc=reshape(permute(v,[1 3 2]),[6.*nc nc]);

% Calculate angles.
norm=sqrt(uZc.*uZc+vZc.*vZc);
AngleCS =  uZc./norm;
AngleSN = -vZc./norm;
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;