Warning, /utils/matlab/cs_grid/uvLatLon2cube.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit a2f8fce7 on 2018-03-10 18:15:33 UTC
f7548778f6 Andr*0001 function [uCs,vCs,errFlag]=uvLatLon2cube(xc,yc,uFld,vFld,xcs,ycs,spv,cosalpha,sinalpha);
a2f8fce7a4 Davi*0002
0003 % [uCs,vCs]=uvLatLon2cube(xc,yc,uFld,vFld,xcs,ycs,spv,cosalpha,sinalpha);
896b20e91e Jean*0004 % put a vector field (uFld,vFld) on to the C-grid, Cubed-sphere grid: uCs,vCs
a2f8fce7a4 Davi*0005 % xc,yc : long,lat position of vector field (uFld,vFld)
896b20e91e Jean*0006 % assume: size(xc) == size(uFld,1) == size(vFld,1)
0007 % and : size(yc) == size(uFld,2) == size(vFld,2)
a2f8fce7a4 Davi*0008 % xcs,ycs : long,lat of the Cell-Center point on CS-grid
896b20e91e Jean*0009 % assume: size(xcs)=size(ycs)=[6*nc nc]=size(uCs)[1:2]=size(vCs)[1:2]
a2f8fce7a4 Davi*0010 % spv : masking value of lat-lon fields
0011 % cosalpha,sinalpha : [optional] cos and sin of angle of cube-sphered grid with
0012 % lat/lon directions, can be obtained from cubeCalcAngle.m
896b20e91e Jean*0013 %
0014 % Written by jmc@ocean.mit.edu, 2005.
a2f8fce7a4 Davi*0015
896b20e91e Jean*0016 if nargin < 7, mask=0 ; else mask=1 ; end
0017 uCs=0; vCs=0; mCsU=0; mCsV=0; errFlag=0;
0018
0019 %Rac='/home/jmc/grid_cs32/';
0020 Rac='grid_files/';
0021
a2f8fce7a4 Davi*0022 nc=size(xcs,2); ncx=6*nc; nPg=ncx*nc; ncp=nc+1;
896b20e91e Jean*0023 nr=size(uFld,3);
0024
0025 fprintf('Entering uvLatLon2cube:');
0026 msk1=ones(size(uFld));
0027 if mask == 1,
0028 fld=ones(size(vFld));
0029 msk1(find(uFld == spv))=0;
0030 fld(find(vFld == spv))=0;
0031 fld=fld-msk1; dd=sum(sum(sum(abs(fld))));
0032 if dd == 0,
0033 fprintf(' mask from uFld & vFld match\n');
0034 else
0035 fprintf(' different uFld & vFld mask in %i pts\n',dd);
0036 errFlag=1; return;
0037 end
0038 uFld(find(msk1 == 0))=0;
0039 vFld(find(msk1 == 0))=0;
0040 else
0041 fprintf(' no mask applied\n');
0042 end
0043
0044 %-- check longitude:
0045 dd=mean(mean(xcs))-mean(xc);
0046 fprintf('mean longitude (data/CS): %5.3f %5.3f and dx= %5.3f \n', ...
0047 mean(xc),mean(mean(xcs)),xc(2)-xc(1));
0048 fprintf(' min max longitude (data) : %5.3f %5.3f \n',xc(1),xc(end));
0049 fprintf(' min max longitude (CSgrid): %5.3f %5.3f \n',min(min(xcs)),max(max(xcs)));
0050 %- add 1 column at the end :
0051 if xc(end) < max(max(xcs)),
0052 dimFld=size(uFld); nx=dimFld(1); nxp=nx+1; dimFld(1)=nxp;
0053 fld=uFld; uFld=zeros(dimFld);
0054 uFld(1:nx,:,:)=fld(1:nx,:,:); uFld(nxp,:,:)=fld(1,:,:);
0055 fld=vFld; vFld=zeros(dimFld);
0056 vFld(1:nx,:,:)=fld(1:nx,:,:); vFld(nxp,:,:)=fld(1,:,:);
0057 fld=msk1; msk1=zeros(dimFld);
0058 msk1(1:nx,:,:)=fld(1:nx,:,:); msk1(nxp,:,:)=fld(1,:,:);
0059 xExt=zeros(1,nxp); xExt(1:nx)=xc; xExt(nxp)=xc(1)+360;
0060 fprintf('Add one column at the end (i=%i): lon= %5.3f\n',dimFld(1),xExt(end));
0061 else
0062 xExt=xc;
0063 end
0064
f7548778f6 Andr*0065 % Read cos and sin of rotation angle if not provided on input
0066 if nargin < 9
0067 namfil=['proj_cs',int2str(nc),'_2uEvN.bin'];
0068 cosalpha=zeros(nPg); sinalpha=zeros(nPg);
0069 fid=fopen([Rac,namfil],'r','b');
0070 cosalpha=fread(fid,nPg,'real*8');
0071 sinalpha=fread(fid,nPg,'real*8');
0072 fclose(fid);
0073 end
0074
0075 uvEN=zeros(nPg, 2);
0076 uvEN(:,1)=cosalpha;
0077 uvEN(:,2)=sinalpha;
896b20e91e Jean*0078
0079 %- go to CS grid, keeping E-W, N-S direction:
0080 x6s=permute(reshape(xcs,[nc 6 nc]),[1 3 2]);
0081 y6s=permute(reshape(ycs,[nc 6 nc]),[1 3 2]);
0082 [uCsE]=int2CS(uFld,xExt,yc,x6s,y6s);
0083 [vCsN]=int2CS(vFld,xExt,yc,x6s,y6s);
0084 [mskC]=int2CS(msk1,xExt,yc,x6s,y6s);
0085 uCsE=uCsE./mskC;
0086 vCsN=vCsN./mskC;
0087 uCsE(find(mskC==0))=0;
0088 vCsN(find(mskC==0))=0;
0089
0090 %- rotate toward CS-grid directions :
0091 uCsE=reshape(permute(uCsE,[1 3 2 4]),nPg,nr);
0092 vCsN=reshape(permute(vCsN,[1 3 2 4]),nPg,nr);
0093
0094 ucs=zeros(nPg,nr);
0095 vcs=zeros(nPg,nr);
0096
0097 for k=1:nr,
0098 ucs(:,k)= uvEN(:,1).*uCsE(:,k)+uvEN(:,2).*vCsN(:,k);
0099 vcs(:,k)=-uvEN(:,2).*uCsE(:,k)+uvEN(:,1).*vCsN(:,k);
0100 end
0101
0102 ucs=reshape(ucs,[ncx nc nr]);
0103 vcs=reshape(vcs,[ncx nc nr]);
0104 mskC=reshape(permute(mskC,[1 3 2 4]),[ncx nc nr]);
0105
0106 %- from A-grid to C-grid :
0107 [u6t]=split_C_cub(ucs);
0108 [v6t]=split_C_cub(vcs);
0109 [m6t]=split_C_cub(mskC);
a2f8fce7a4 Davi*0110 if nr == 1
0111 u6t = reshape(u6t,ncp,ncp,nr,6);
0112 v6t = reshape(v6t,ncp,ncp,nr,6);
0113 m6t = reshape(m6t,ncp,ncp,nr,6);
0114 end
896b20e91e Jean*0115 for n=1:2:6,
0116 uloc=u6t(1,:,:,n);
0117 u6t(1,:,:,n)=v6t(1,:,:,n);
0118 v6t(1,:,:,n)=uloc;
0119 end
0120 for n=2:2:6,
0121 uloc=u6t(:,1,:,n);
0122 u6t(:,1,:,n)=v6t(:,1,:,n);
0123 v6t(:,1,:,n)=uloc;
0124 end
0125
0126 uCs=zeros(nc,nc,nr,6);
0127 vCs=zeros(nc,nc,nr,6);
0128 mCsU=zeros(nc,nc,nr,6);
0129 mCsV=zeros(nc,nc,nr,6);
0130 uCs(:,:,:,:)=u6t([1:nc],[2:ncp],:,:)+u6t([2:ncp],[2:ncp],:,:);
0131 mCsU(:,:,:,:)=m6t([1:nc],[2:ncp],:,:)+m6t([2:ncp],[2:ncp],:,:);
0132 vCs(:,:,:,:)=v6t([2:ncp],[1:nc],:,:)+v6t([2:ncp],[2:ncp],:,:);
0133 mCsV(:,:,:,:)=m6t([2:ncp],[1:nc],:,:)+m6t([2:ncp],[2:ncp],:,:);
0134 uCs=reshape(permute(uCs,[1 4 2 3]),[ncx nc nr])/2;
0135 vCs=reshape(permute(vCs,[1 4 2 3]),[ncx nc nr])/2;
0136 mCsU=reshape(permute(mCsU,[1 4 2 3]),[ncx nc nr])/2;
0137 mCsV=reshape(permute(mCsV,[1 4 2 3]),[ncx nc nr])/2;
0138
0139 if mask == 1,
0140 uCs(find(mCsU==0.))=spv;
0141 vCs(find(mCsV==0.))=spv;
0142 end
0143
0144 return
0145
0146 function Qcs=int2CS(Qini, x0,y0,xc,yc)
0147
0148 Qcs=zeros([size(xc) size(Qini,3)]);
0149 for m=1:size(Qini,3)
0150 for f=1:size(xc,3);
0151 Qcs(:,:,f,m)=interp2(y0,x0,Qini(:,:,m),yc(:,:,f),xc(:,:,f));
0152 end
0153 end
0154
0155 return
