Warning, /utils/matlab/cs_grid/calc2ndmom.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
70117b93f1 Andr*0001 function [UVtot,UVtrans,U2tot,U2trans,V2tot,V2trans,errFlag]=calc2ndmom(u,v,usq,vsq,uv,cosalpha,sinalpha);
896b20e91e Jean*0002 % [UVtot,UVtrans,U2tot,U2trans,V2tot,V2trans,errFlag]=calc2ndmom(u,v,usq,vsq,uv);
0003 % Compute UV, U2 and V2 (at cubed-sphere A grid points) using cubed-sphere input
0004 % u = time averaged u-wind on cubed sphere grid at u-points
0005 % v = time averaged v-wind on cubed sphere grid at v-points
0006 % usq = time averaged square of u-wind on cubed sphere grid at u-points
0007 % vsq = time averaged square of v-wind on cubed sphere grid at v-points
0008 % uv = time averaged product of u-wind and v-wind on cubed sphere grid at mass points
0009 %
0010 % Written by molod@ocean.mit.edu, 2005.
0011 fprintf('Entering calc2ndmom: \n');
0012 errFlag=0;
0013 %
0014 nc=size(u,2); ncx=6*nc; nPg=ncx*nc;
0015 nz=size(u,3);
0016 %
0017 %Check size of input arrays. if read by rdmnc, remove column and row
0018 % if read by rdmds, leave alone
0019 % if wrong size altogether, get out with error
0020 %
0021 NN = size(u);
0022 if NN(1) == 6*nc+1,
0023 fprintf('Need to remove a column \n');
05128db2bc Davi*0024 U = u(1:6*nc,:);
0025 U2 = usq(1:6*nc,:);
896b20e91e Jean*0026 elseif NN(1) == 6*nc,
0027 fprintf('No U reshaping needed \n');
05128db2bc Davi*0028 U = u;
0029 U2 = usq;
896b20e91e Jean*0030 else
0031 fprintf('Error in U-point CS-dim: %i %i %i \n',NN);
0032 return
0033 end
0034 %
0035 NN = size(v);
0036 if NN(2) == nc+1,
0037 fprintf('Need to remove a row ');
05128db2bc Davi*0038 V = v(:,1:nc);
0039 V2 = vsq(:,1:nc);
896b20e91e Jean*0040 elseif NN(2) == nc,
0041 fprintf('No V reshaping needed \n');
05128db2bc Davi*0042 V = v;
0043 V2 = vsq;
896b20e91e Jean*0044 else
0045 fprintf('Error in V-point CS-dim: %i %i %i \n',NN);
0046 return
0047 end
0048 %
0049 %Get ready to do the exchange
0050 %
05128db2bc Davi*0051 U=reshape(U,[nc 6 nc nz]);
0052 U2=reshape(U2,[nc 6 nc nz]);
0053 V=reshape(V,[nc 6 nc nz]);
0054 V2=reshape(V2,[nc 6 nc nz]);
896b20e91e Jean*0055 %
0056 uu=zeros(nc+1,6,nc,nz);
0057 uu2=zeros(nc+1,6,nc,nz);
0058 vv=zeros(nc,6,nc+1,nz);
0059 vv2=zeros(nc,6,nc+1,nz);
0060 %
0061 for ifc=1:6;
0062 uu(1:nc,ifc,:,:)=U(:,ifc,:,:);
0063 uu2(1:nc,ifc,:,:)=U2(:,ifc,:,:);
0064 vv(:,ifc,1:nc,:)=V(:,ifc,:,:);
0065 vv2(:,ifc,1:nc,:)=V2(:,ifc,:,:);
0066 end
0067 %
0068 % do the exchange
0069 for k=1:nz;
0070 uu(nc+1,1:2:6,:,k)=uu(1,2:2:6,:,k);
0071 uu(nc+1,2:2:6,:,k)=vv(nc:-1:1,[4:2:6 2:2:3],1,k)';
0072 uu2(nc+1,1:2:6,:,k)=uu2(1,2:2:6,:,k);
0073 uu2(nc+1,2:2:6,:,k)=vv2(nc:-1:1,[4:2:6 2:2:3],1,k)';
0074 vv(:,2:2:6,nc+1,k)=vv(:,[3:2:6 1],1,k);
0075 vv(:,1:2:6,nc+1,k)=squeeze(uu(1,[3:2:6 1],nc:-1:1,k))';
0076 vv2(:,2:2:6,nc+1,k)=vv2(:,[3:2:6 1],1,k);
0077 vv2(:,1:2:6,nc+1,k)=squeeze(uu2(1,[3:2:6 1],nc:-1:1,k))';
0078 end
0079 %
0080 % Interp to mass points
0081 ui=(uu(1:nc,:,:,:)+uu(2:nc+1,:,:,:))/2;
0082 ub2i=(uu(1:nc,:,:,:).*uu(1:nc,:,:,:)+uu(2:nc+1,:,:,:).*uu(2:nc+1,:,:,:))/2;
0083 u2i=(uu2(1:nc,:,:,:)+uu2(2:nc+1,:,:,:))/2;
0084 vj=(vv(:,:,1:nc,:)+vv(:,:,2:nc+1,:))/2;
0085 vb2j=(vv(:,:,1:nc,:).*vv(:,:,1:nc,:)+vv(:,:,2:nc+1,:).*vv(:,:,2:nc+1,:))/2;
0086 v2j=(vv2(:,:,1:nc,:)+vv2(:,:,2:nc+1,:))/2;
0087 %
0088 ui=reshape(ui,[nPg nz]);
0089 u2i=reshape(u2i,[nPg nz]);
6f3388b7c3 Andr*0090 ub2i=reshape(ub2i,[nPg nz]);
896b20e91e Jean*0091 vj=reshape(vj,[nPg nz]);
0092 v2j=reshape(v2j,[nPg nz]);
6f3388b7c3 Andr*0093 vb2j=reshape(vb2j,[nPg nz]);
896b20e91e Jean*0094 UV=reshape(uv,[nPg nz]);
0095 %
70117b93f1 Andr*0096 % Read cos and sin of rotation angle if not provided on input
0097 if nargin == 5
0098 Rac='/u/u2/jmc/';
0099 namfil=['proj_cs',int2str(nc),'_2uEvN.bin'];
0100 cosalpha=zeros(nPg); sinalpha=zeros(nPg);
0101 fid=fopen([Rac,namfil],'r','b');
0102 cosalpha=fread(fid,nPg,'real*8');
0103 sinalpha=fread(fid,nPg,'real*8');
0104 fclose(fid);
0105 end
896b20e91e Jean*0106 %
0107 UVtrans=zeros(nPg,nz); UVtot=zeros(nPg,nz);
0108 U2trans=zeros(nPg,nz); U2tot=zeros(nPg,nz);
0109 V2trans=zeros(nPg,nz); V2tot=zeros(nPg,nz);
0110 %
0111 for k=1:nz;
0112 UVtrans(:,k)=( (u2i(:,k)-ub2i(:,k)) - (v2j(:,k)-vb2j(:,k))) .* cosalpha .* sinalpha ...
0113 + ( UV(:,k) - ui(:,k).*vj(:,k) ) .* (cosalpha.*cosalpha - sinalpha.*sinalpha);
0114 UVtot(:,k)=UVtrans(:,k) + ( ui(:,k).*ui(:,k) - vj(:,k).*vj(:,k) ) .* cosalpha .* sinalpha ...
0115 + ui(:,k).*vj(:,k).*(cosalpha.*cosalpha - sinalpha.*sinalpha);
0116 U2tot(:,k) = u2i(:,k) .* cosalpha.*cosalpha + v2j(:,k).*sinalpha.*sinalpha ...
0117 - 2.* UV(:,k) .*cosalpha .* sinalpha;
0118 U2trans(:,k) = U2tot(:,k) - ...
0119 ( ui(:,k).*ui(:,k).*cosalpha.*cosalpha + vj(:,k).*vj(:,k).*sinalpha.*sinalpha ...
0120 - 2.* ui(:,k).*vj(:,k) .*cosalpha .* sinalpha);
0121 V2tot(:,k) = u2i(:,k) .* sinalpha.*sinalpha + v2j(:,k).*cosalpha.*cosalpha ...
0122 + 2.* UV(:,k) .*cosalpha .* sinalpha;
0123 V2trans(:,k) = V2tot(:,k) - ...
0124 ( ui(:,k).*ui(:,k).*sinalpha.*sinalpha + vj(:,k).*vj(:,k).*cosalpha.*cosalpha ...
0125 + 2.* ui(:,k).*vj(:,k) .*cosalpha .* sinalpha);
0126 end
0127 %
0128 UVtot=reshape(UVtot,[6*nc nc nz]);
0129 UVtrans=reshape(UVtrans,[6*nc nc nz]);
0130 U2tot=reshape(U2tot,[6*nc nc nz]);
0131 U2trans=reshape(U2trans,[6*nc nc nz]);
0132 V2tot=reshape(V2tot,[6*nc nc nz]);
0133 V2trans=reshape(V2trans,[6*nc nc nz]);
0134 %
0135 clear sinalpha cosalpha
