| 1 |
function [uCs,vCs,errFlag]=uvLatLon2cube(xc,yc,uFld,vFld,xcs,ycs,spv); |
| 2 |
% [uCs,vCs]=uvLatLon2cube(xc,yc,uFld,vFld,xcs,ycs,spv); |
| 3 |
% put a vector field (uFld,vFld) on to the C-grid, Cubed-sphere grid: uCs,vCs |
| 4 |
% xc,yc = long,lat position of vector field (uFld,vFld) |
| 5 |
% assume: size(xc) == size(uFld,1) == size(vFld,1) |
| 6 |
% and : size(yc) == size(uFld,2) == size(vFld,2) |
| 7 |
% xcs,ycs = long,lat of the Cell-Center point on CS-grid |
| 8 |
% assume: size(xcs)=size(ycs)=[6*nc nc]=size(uCs)[1:2]=size(vCs)[1:2] |
| 9 |
% $Header: /u/gcmpack/MITgcm/utils/cs_grid/uvLatLon2cube.m,v 1.2 2005/08/16 20:34:02 molod Exp $ |
| 10 |
|
| 11 |
if nargin < 7, mask=0 ; else mask=1 ; end |
| 12 |
uCs=0; vCs=0; mCsU=0; mCsV=0; errFlag=0; |
| 13 |
|
| 14 |
%Rac='/home/jmc/grid_cs32/'; |
| 15 |
Rac='grid_files/'; |
| 16 |
|
| 17 |
nc=size(xcs,2); ncx=6*nc; nPg=ncx*nc; |
| 18 |
nr=size(uFld,3); |
| 19 |
|
| 20 |
fprintf('Entering uvLatLon2cube:'); |
| 21 |
msk1=ones(size(uFld)); |
| 22 |
if mask == 1, |
| 23 |
fld=ones(size(vFld)); |
| 24 |
msk1(find(uFld == spv))=0; |
| 25 |
fld(find(vFld == spv))=0; |
| 26 |
fld=fld-msk1; dd=sum(sum(sum(abs(fld)))); |
| 27 |
if dd == 0, |
| 28 |
fprintf(' mask from uFld & vFld match\n'); |
| 29 |
else |
| 30 |
fprintf(' different uFld & vFld mask in %i pts\n',dd); |
| 31 |
errFlag=1; return; |
| 32 |
end |
| 33 |
uFld(find(msk1 == 0))=0; |
| 34 |
vFld(find(msk1 == 0))=0; |
| 35 |
else |
| 36 |
fprintf(' no mask applied\n'); |
| 37 |
end |
| 38 |
|
| 39 |
%-- check longitude: |
| 40 |
dd=mean(mean(xcs))-mean(xc); |
| 41 |
fprintf('mean longitude (data/CS): %5.3f %5.3f and dx= %5.3f \n', ... |
| 42 |
mean(xc),mean(mean(xcs)),xc(2)-xc(1)); |
| 43 |
fprintf(' min max longitude (data) : %5.3f %5.3f \n',xc(1),xc(end)); |
| 44 |
fprintf(' min max longitude (CSgrid): %5.3f %5.3f \n',min(min(xcs)),max(max(xcs))); |
| 45 |
%- add 1 column at the end : |
| 46 |
if xc(end) < max(max(xcs)), |
| 47 |
dimFld=size(uFld); nx=dimFld(1); nxp=nx+1; dimFld(1)=nxp; |
| 48 |
fld=uFld; uFld=zeros(dimFld); |
| 49 |
uFld(1:nx,:,:)=fld(1:nx,:,:); uFld(nxp,:,:)=fld(1,:,:); |
| 50 |
fld=vFld; vFld=zeros(dimFld); |
| 51 |
vFld(1:nx,:,:)=fld(1:nx,:,:); vFld(nxp,:,:)=fld(1,:,:); |
| 52 |
fld=msk1; msk1=zeros(dimFld); |
| 53 |
msk1(1:nx,:,:)=fld(1:nx,:,:); msk1(nxp,:,:)=fld(1,:,:); |
| 54 |
xExt=zeros(1,nxp); xExt(1:nx)=xc; xExt(nxp)=xc(1)+360; |
| 55 |
fprintf('Add one column at the end (i=%i): lon= %5.3f\n',dimFld(1),xExt(end)); |
| 56 |
else |
| 57 |
xExt=xc; |
| 58 |
end |
| 59 |
|
| 60 |
namfil=['proj_cs',int2str(nc),'_2uEvN.bin']; |
| 61 |
fid=fopen([Rac,namfil],'r','b'); uvEN=fread(fid,nPg*2,'real*8'); fclose(fid); |
| 62 |
uvEN=reshape(uvEN,[nPg 2]); |
| 63 |
|
| 64 |
%- go to CS grid, keeping E-W, N-S direction: |
| 65 |
x6s=permute(reshape(xcs,[nc 6 nc]),[1 3 2]); |
| 66 |
y6s=permute(reshape(ycs,[nc 6 nc]),[1 3 2]); |
| 67 |
[uCsE]=int2CS(uFld,xExt,yc,x6s,y6s); |
| 68 |
[vCsN]=int2CS(vFld,xExt,yc,x6s,y6s); |
| 69 |
[mskC]=int2CS(msk1,xExt,yc,x6s,y6s); |
| 70 |
uCsE=uCsE./mskC; |
| 71 |
vCsN=vCsN./mskC; |
| 72 |
uCsE(find(mskC==0))=0; |
| 73 |
vCsN(find(mskC==0))=0; |
| 74 |
|
| 75 |
%- rotate toward CS-grid directions : |
| 76 |
uCsE=reshape(permute(uCsE,[1 3 2 4]),nPg,nr); |
| 77 |
vCsN=reshape(permute(vCsN,[1 3 2 4]),nPg,nr); |
| 78 |
|
| 79 |
ucs=zeros(nPg,nr); |
| 80 |
vcs=zeros(nPg,nr); |
| 81 |
|
| 82 |
for k=1:nr, |
| 83 |
ucs(:,k)= uvEN(:,1).*uCsE(:,k)+uvEN(:,2).*vCsN(:,k); |
| 84 |
vcs(:,k)=-uvEN(:,2).*uCsE(:,k)+uvEN(:,1).*vCsN(:,k); |
| 85 |
end |
| 86 |
|
| 87 |
ucs=reshape(ucs,[ncx nc nr]); |
| 88 |
vcs=reshape(vcs,[ncx nc nr]); |
| 89 |
mskC=reshape(permute(mskC,[1 3 2 4]),[ncx nc nr]); |
| 90 |
|
| 91 |
%- from A-grid to C-grid : |
| 92 |
[u6t]=split_C_cub(ucs); |
| 93 |
[v6t]=split_C_cub(vcs); |
| 94 |
[m6t]=split_C_cub(mskC); |
| 95 |
for n=1:2:6, |
| 96 |
uloc=u6t(1,:,:,n); |
| 97 |
u6t(1,:,:,n)=v6t(1,:,:,n); |
| 98 |
v6t(1,:,:,n)=uloc; |
| 99 |
end |
| 100 |
for n=2:2:6, |
| 101 |
uloc=u6t(:,1,:,n); |
| 102 |
u6t(:,1,:,n)=v6t(:,1,:,n); |
| 103 |
v6t(:,1,:,n)=uloc; |
| 104 |
end |
| 105 |
%size(u6t) |
| 106 |
ncp=nc+1; |
| 107 |
|
| 108 |
uCs=zeros(nc,nc,nr,6); |
| 109 |
vCs=zeros(nc,nc,nr,6); |
| 110 |
mCsU=zeros(nc,nc,nr,6); |
| 111 |
mCsV=zeros(nc,nc,nr,6); |
| 112 |
uCs(:,:,:,:)=u6t([1:nc],[2:ncp],:,:)+u6t([2:ncp],[2:ncp],:,:); |
| 113 |
mCsU(:,:,:,:)=m6t([1:nc],[2:ncp],:,:)+m6t([2:ncp],[2:ncp],:,:); |
| 114 |
vCs(:,:,:,:)=v6t([2:ncp],[1:nc],:,:)+v6t([2:ncp],[2:ncp],:,:); |
| 115 |
mCsV(:,:,:,:)=m6t([2:ncp],[1:nc],:,:)+m6t([2:ncp],[2:ncp],:,:); |
| 116 |
uCs=reshape(permute(uCs,[1 4 2 3]),[ncx nc nr])/2; |
| 117 |
vCs=reshape(permute(vCs,[1 4 2 3]),[ncx nc nr])/2; |
| 118 |
mCsU=reshape(permute(mCsU,[1 4 2 3]),[ncx nc nr])/2; |
| 119 |
mCsV=reshape(permute(mCsV,[1 4 2 3]),[ncx nc nr])/2; |
| 120 |
|
| 121 |
if mask == 1, |
| 122 |
uCs(find(mCsU==0.))=spv; |
| 123 |
vCs(find(mCsV==0.))=spv; |
| 124 |
end |
| 125 |
|
| 126 |
return |
| 127 |
|
| 128 |
function Qcs=int2CS(Qini, x0,y0,xc,yc) |
| 129 |
|
| 130 |
Qcs=zeros([size(xc) size(Qini,3)]); |
| 131 |
for m=1:size(Qini,3) |
| 132 |
for f=1:size(xc,3); |
| 133 |
Qcs(:,:,f,m)=interp2(y0,x0,Qini(:,:,m),yc(:,:,f),xc(:,:,f)); |
| 134 |
end |
| 135 |
end |
| 136 |
|
| 137 |
return |
Parent Directory
| 
Revision Log
| 
Revision Graph