| 1 |
function [UVtot,UVtrans,U2tot,U2trans,V2tot,V2trans,errFlag]=calc2ndmom(u,v,usq,vsq,uv); |
| 2 |
% [UVtot,UVtrans,U2tot,U2trans,V2tot,V2trans,errFlag]=calc2ndmom(u,v,usq,vsq,uv); |
| 3 |
% Compute UV, U2 and V2 (at cubed-sphere A grid points) using cubed-sphere input |
| 4 |
% u = time averaged u-wind on cubed sphere grid at u-points |
| 5 |
% v = time averaged v-wind on cubed sphere grid at v-points |
| 6 |
% usq = time averaged square of u-wind on cubed sphere grid at u-points |
| 7 |
% vsq = time averaged square of v-wind on cubed sphere grid at v-points |
| 8 |
% uv = time averaged product of u-wind and v-wind on cubed sphere grid at mass points |
| 9 |
% |
| 10 |
fprintf('Entering calc2ndmom: \n'); |
| 11 |
errFlag=0; |
| 12 |
% |
| 13 |
nc=size(u,2); ncx=6*nc; nPg=ncx*nc; |
| 14 |
nz=size(u,3); |
| 15 |
% |
| 16 |
%Check size of input arrays. if read by rdmnc, remove column and row |
| 17 |
% if read by rdmds, leave alone |
| 18 |
% if wrong size altogether, get out with error |
| 19 |
% |
| 20 |
NN = size(u); |
| 21 |
if NN(1) == 6*nc+1, |
| 22 |
fprintf('Need to remove a column \n'); |
| 23 |
utmp = u(1:6*nc,:); |
| 24 |
u2tmp = usq(1:6*nc,:); |
| 25 |
elseif NN(1) == 6*nc, |
| 26 |
fprintf('No U reshaping needed \n'); |
| 27 |
utmp = u; |
| 28 |
u2tmp = usq; |
| 29 |
else |
| 30 |
fprintf('Error in U-point CS-dim: %i %i %i \n',NN); |
| 31 |
return |
| 32 |
end |
| 33 |
% |
| 34 |
NN = size(v); |
| 35 |
if NN(2) == nc+1, |
| 36 |
fprintf('Need to remove a row '); |
| 37 |
vtmp = v(:,1:nc); |
| 38 |
v2tmp = vsq(:,1:nc); |
| 39 |
elseif NN(2) == nc, |
| 40 |
fprintf('No V reshaping needed \n'); |
| 41 |
vtmp = v; |
| 42 |
v2tmp = vsq; |
| 43 |
else |
| 44 |
fprintf('Error in V-point CS-dim: %i %i %i \n',NN); |
| 45 |
return |
| 46 |
end |
| 47 |
% |
| 48 |
%Get ready to do the exchange |
| 49 |
% |
| 50 |
U=reshape(utmp,[nc 6 nc nz]); |
| 51 |
U2=reshape(u2tmp,[nc 6 nc nz]); |
| 52 |
V=reshape(vtmp,[nc 6 nc nz]); |
| 53 |
V2=reshape(v2tmp,[nc 6 nc nz]); |
| 54 |
% |
| 55 |
uu=zeros(nc+1,6,nc,nz); |
| 56 |
uu2=zeros(nc+1,6,nc,nz); |
| 57 |
vv=zeros(nc,6,nc+1,nz); |
| 58 |
vv2=zeros(nc,6,nc+1,nz); |
| 59 |
% |
| 60 |
for ifc=1:6; |
| 61 |
uu(1:nc,ifc,:,:)=U(:,ifc,:,:); |
| 62 |
uu2(1:nc,ifc,:,:)=U2(:,ifc,:,:); |
| 63 |
vv(:,ifc,1:nc,:)=V(:,ifc,:,:); |
| 64 |
vv2(:,ifc,1:nc,:)=V2(:,ifc,:,:); |
| 65 |
end |
| 66 |
% |
| 67 |
% do the exchange |
| 68 |
for k=1:nz; |
| 69 |
uu(nc+1,1:2:6,:,k)=uu(1,2:2:6,:,k); |
| 70 |
uu(nc+1,2:2:6,:,k)=vv(nc:-1:1,[4:2:6 2:2:3],1,k)'; |
| 71 |
uu2(nc+1,1:2:6,:,k)=uu2(1,2:2:6,:,k); |
| 72 |
uu2(nc+1,2:2:6,:,k)=vv2(nc:-1:1,[4:2:6 2:2:3],1,k)'; |
| 73 |
vv(:,2:2:6,nc+1,k)=vv(:,[3:2:6 1],1,k); |
| 74 |
vv(:,1:2:6,nc+1,k)=squeeze(uu(1,[3:2:6 1],nc:-1:1,k))'; |
| 75 |
vv2(:,2:2:6,nc+1,k)=vv2(:,[3:2:6 1],1,k); |
| 76 |
vv2(:,1:2:6,nc+1,k)=squeeze(uu2(1,[3:2:6 1],nc:-1:1,k))'; |
| 77 |
end |
| 78 |
% |
| 79 |
% Interp to mass points |
| 80 |
ui=(uu(1:nc,:,:,:)+uu(2:nc+1,:,:,:))/2; |
| 81 |
ub2i=(uu(1:nc,:,:,:).*uu(1:nc,:,:,:)+uu(2:nc+1,:,:,:).*uu(2:nc+1,:,:,:))/2; |
| 82 |
u2i=(uu2(1:nc,:,:,:)+uu2(2:nc+1,:,:,:))/2; |
| 83 |
vj=(vv(:,:,1:nc,:)+vv(:,:,2:nc+1,:))/2; |
| 84 |
vb2j=(vv(:,:,1:nc,:).*vv(:,:,1:nc,:)+vv(:,:,2:nc+1,:).*vv(:,:,2:nc+1,:))/2; |
| 85 |
v2j=(vv2(:,:,1:nc,:)+vv2(:,:,2:nc+1,:))/2; |
| 86 |
% |
| 87 |
ui=reshape(ui,[nPg nz]); |
| 88 |
u2i=reshape(u2i,[nPg nz]); |
| 89 |
ub2i=reshape(u2i,[nPg nz]); |
| 90 |
vj=reshape(vj,[nPg nz]); |
| 91 |
v2j=reshape(v2j,[nPg nz]); |
| 92 |
vb2j=reshape(v2j,[nPg nz]); |
| 93 |
UV=reshape(uv,[nPg nz]); |
| 94 |
% |
| 95 |
% Read cos and sin of rotation angle |
| 96 |
Rac='/u/u2/jmc/'; |
| 97 |
namfil=['proj_cs',int2str(nc),'_2uEvN.bin']; |
| 98 |
cosalpha=zeros(nPg); sinalpha=zeros(nPg); |
| 99 |
fid=fopen([Rac,namfil],'r','b'); |
| 100 |
cosalpha=fread(fid,nPg,'real*8'); |
| 101 |
sinalpha=fread(fid,nPg,'real*8'); |
| 102 |
fclose(fid); |
| 103 |
% |
| 104 |
UVtrans=zeros(nPg,nz); UVtot=zeros(nPg,nz); |
| 105 |
U2trans=zeros(nPg,nz); U2tot=zeros(nPg,nz); |
| 106 |
V2trans=zeros(nPg,nz); V2tot=zeros(nPg,nz); |
| 107 |
% |
| 108 |
for k=1:nz; |
| 109 |
UVtrans(:,k)=( (u2i(:,k)-ub2i(:,k)) - (v2j(:,k)-vb2j(:,k))) .* cosalpha .* sinalpha ... |
| 110 |
+ ( UV(:,k) - ui(:,k).*vj(:,k) ) .* (cosalpha.*cosalpha - sinalpha.*sinalpha); |
| 111 |
UVtot(:,k)=UVtrans(:,k) + ( ui(:,k).*ui(:,k) - vj(:,k).*vj(:,k) ) .* cosalpha .* sinalpha ... |
| 112 |
+ ui(:,k).*vj(:,k).*(cosalpha.*cosalpha - sinalpha.*sinalpha); |
| 113 |
% UVtot(:,k) = (u2i(:,k)-v2j(:,k)).* sinalpha.*cosalpha + ... |
| 114 |
% UV(:,k).*(cosalpha.*cosalpha - sinalpha.*sinalpha); |
| 115 |
% UVtrans(:,k) = UVtot(:,k) - ... |
| 116 |
% ( (ui(:,k).*ui(:,k)-vj(:,k).*vj(:,k)).* sinalpha.*cosalpha + ... |
| 117 |
% ui(:,k).*vj(:,k).*(cosalpha.*cosalpha - sinalpha.*sinalpha) ); |
| 118 |
U2tot(:,k) = u2i(:,k) .* cosalpha.*cosalpha + v2j(:,k).*sinalpha.*sinalpha ... |
| 119 |
- 2.* UV(:,k) .*cosalpha .* sinalpha; |
| 120 |
U2trans(:,k) = U2tot(:,k) - ... |
| 121 |
( ui(:,k).*ui(:,k).*cosalpha.*cosalpha + vj(:,k).*vj(:,k).*sinalpha.*sinalpha ... |
| 122 |
- 2.* ui(:,k).*vj(:,k) .*cosalpha .* sinalpha); |
| 123 |
V2tot(:,k) = u2i(:,k) .* sinalpha.*sinalpha + v2j(:,k).*cosalpha.*cosalpha ... |
| 124 |
+ 2.* UV(:,k) .*cosalpha .* sinalpha; |
| 125 |
V2trans(:,k) = V2tot(:,k) - ... |
| 126 |
( ui(:,k).*ui(:,k).*sinalpha.*sinalpha + vj(:,k).*vj(:,k).*cosalpha.*cosalpha ... |
| 127 |
+ 2.* ui(:,k).*vj(:,k) .*cosalpha .* sinalpha); |
| 128 |
end |
| 129 |
% |
| 130 |
UVtot=reshape(UVtot,[6*nc nc nz]); |
| 131 |
UVtrans=reshape(UVtrans,[6*nc nc nz]); |
| 132 |
U2tot=reshape(U2tot,[6*nc nc nz]); |
| 133 |
U2trans=reshape(U2trans,[6*nc nc nz]); |
| 134 |
V2tot=reshape(V2tot,[6*nc nc nz]); |
| 135 |
V2trans=reshape(V2trans,[6*nc nc nz]); |
| 136 |
% |
| 137 |
clear sinalpha cosalpha |
Parent Directory
| 
Revision Log
| 
Revision Graph