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; |