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

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

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


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

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

  ViewVC Help
Powered by ViewVC 1.1.22