/[MITgcm]/MITgcm/utils/cs_grid/calc2ndmom.m
ViewVC logotype

Contents of /MITgcm/utils/cs_grid/calc2ndmom.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (show annotations) (download)
Thu Sep 15 16:51:13 2005 UTC (17 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
moved from utils/cs_grid to utils/matlab/cs_grid.

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

  ViewVC Help
Powered by ViewVC 1.1.22