Warning, /utils/matlab/cs_grid/mk_isoLat_bkl.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 4ec37fd8 on 2023-10-03 22:00:32 UTC
4ec37fd829 Jean*0001 % script to generate broken lines that folows the cube-sphere grid
0002 % and stay close as possible to latitude-target "yLat":
0003 %- method:
0004 % Tag grid-cell center pts (fct of where lat is within yLat separation)
0005 % and set bk-line where 2 adjacent tagC match a crossing of yLat boundary
0006 %-----------------------
0007 %- load definition of the grid (needs to be done at the 1rst call):
0008
0009 if size(who('krd'),1) > 0,
0010 fprintf('krd is defined and = %i \n',krd);
0011 else
0012 fprintf('krd undefined ; set to 1 \n'); krd=1 ;
0013 end
0014 if krd==1,
0015 more off ;
0016 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0017 % or ncdf=0 to load MDS (binary) grid-files :
0018 ncdf=0;
0019 gDir='grid_files/';
0020 G=load_grid(gDir,10+ncdf);
0021 xcs=G.xC; ycs=G.yC; xcg=G.xG; ycg=G.yG;
0022 %- 1rst try: skip writing to file
0023 kwr=0; kplot=2;
0024 else
0025 %- write/save to file:
0026 kwr=1; kplot=0;
0027 end
0028 %- output arrays:
0029 yLat=-87:3:87;
0030 %yLat=-88:2:88;
0031 ydim=length(yLat);
0032
0033 %------------
0034 n1h=size(xcs,1); n2h=size(xcs,2);
0035 if n1h == 6*n2h, nc=n2h;
0036 elseif n1h*6 == n2h, nc=n1h;
0037 else
0038 error([' grid var size: ',int2str(n1h),' x ',int2str(n2h),' does not fit regular cube !']);
0039 end
0040 nPg=nc*nc*6; ncp=nc+1; nPp2=nPg+2;
0041
0042 %- storage arrays:
0043 savNpts=zeros(1,ydim);
0044 savFlag=zeros(6*ncp,ydim);
0045 savIJuv=zeros(6*ncp,ydim);
0046
0047 %-------------------------------------------------------------------------------
0048 %- A) first step: Tag each grid-cell center point
0049 % tagC(i,j) = jl if yLat(jl) <= yC(i,j) < yLat(jl+1)
0050
0051 tagC=-ones(n1h,n2h);
0052 yFull=[-90 yLat +90];
0053 for jl=0:ydim,
0054 %- there is a shift of 1 between yFull and yLat so that yLat(jl) = yFull(jl+1)
0055 var=(ycs-yFull(jl+1)).*(yFull(jl+2)-ycs);
0056 tagC(find(var > 0))=jl;
0057 tagC(find(ycs == yFull(jl+1)))=jl;
0058 end
0059 %tagC(find(ycs < yLat(1)))=0;
0060 %tagC(find(ycs >= yLat(end)))=ydim;
0061
0062 %- B) second step: split to 6 faces and check how big the jump is from 2 neighbor cells
0063 tagC6f = split_C_cub(tagC);
0064
0065 MxJump=0;
0066 nMxdbg=10; %- Max number of debug-print per face
0067 for n=1:6,
0068 var=tagC6f(1:nc,2:ncp,n)-tagC6f(2:ncp,2:ncp,n); var=abs(var);
0069 jumpUpt=max(var(:));
0070 if jumpUpt > 1,
0071 [I J]=find(var > 1);
0072 nbJu=length(I);
0073 fprintf('-- face # %i , %i points with a tagC jump across U.pt:\n',n,nbJu);
0074 for il=1:min(nbJu,nMxdbg),
0075 i=I(il); j=J(il)+1;
0076 fprintf(' -i, j = %i, %i ; x,y= %8.3f, %8.3f ; tagC= %i ( %i )\n', ...
0077 i-1,j-1,xx1(i,j,n),yy1(i,j,n),tagC6f(i,j,n),var(i,j-1))
0078 i=I(il)+1; j=J(il)+1;
0079 fprintf(' +i, j = %i, %i ; x,y= %8.3f, %8.3f ; tagC= %i ( %i )\n', ...
0080 i-1,j-1,xx1(i,j,n),yy1(i,j,n),tagC6f(i,j,n),var(i-1,j-1))
0081 end
0082 end
0083 var=tagC6f(2:ncp,1:nc,n)-tagC6f(2:ncp,2:ncp,n); var=abs(var);
0084 jumpVpt=max(var(:));
0085 if jumpVpt > 1,
0086 [I J]=find(var > 1);
0087 nbJv=length(I);
0088 fprintf('-- face # %i , %i points with a tagC jump across V.pt:\n',n,nbJv);
0089 for il=1:min(nbJv,nMxdbg),
0090 i=I(il)+1; j=J(il);
0091 fprintf(' i, -j = %i, %i ; x,y= %8.3f, %8.3f ; tagC= %i ( %i )\n', ...
0092 i-1,j-1,xx1(i,j,n),yy1(i,j,n),tagC6f(i,j,n),var(i-1,j))
0093 i=I(il)+1; j=J(il)+1;
0094 fprintf(' i, +j = %i, %i ; x,y= %8.3f, %8.3f ; tagC= %i ( %i )\n', ...
0095 i-1,j-1,xx1(i,j,n),yy1(i,j,n),tagC6f(i,j,n),var(i-1,j-1))
0096 end
0097 end
0098 MxJump=max( MxJump, jumpUpt );
0099 MxJump=max( MxJump, jumpVpt );
0100 end
0101 fprintf(' max tagC Jump between 2 neighbor = %i\n',MxJump);
0102
0103 if MxJump == 1,
0104 %===============================================================================
0105 % Case with no Jump > 1 in adjacent tagC:
0106 %-------------------------------------------------------------------------------
0107 %- C.1) third step: Tag U & V point all in one pass:
0108 % tagU(i,j) = +jl if tagC6f(i,j)= jl & tagC6f(i-1,j)= jl-1
0109 % tagU(i,j) = -jl if tagC6f(i,j)= jl-1 & tagC6f(i-1,j)= jl
0110 %
0111 % tagV(i,j) = +jl if tagC6f(i,j)= jl & tagC6f(i,j-1)= jl-1
0112 % tagV(i,j) = -jl if tagC6f(i,j)= jl-1 & tagC6f(i,j-1)= jl
0113 % Note: this way we recover tagU(i,j) == +/- jj if yLat(jl) == Latitude-of-U-pt(i,j)
0114
0115 tagU=zeros(nc,nc,6); tagV=tagU;
0116 for jl=1:ydim,
0117 var=abs(tagC6f(1:nc,2:ncp,:)-jl+1)+abs(tagC6f(2:ncp,2:ncp,:)-jl);
0118 tagU(var==0)=jl;
0119 var=abs(tagC6f(1:nc,2:ncp,:)-jl)+abs(tagC6f(2:ncp,2:ncp,:)-jl+1);
0120 tagU(var==0)=-jl;
0121 var=abs(tagC6f(2:ncp,1:nc,:)-jl+1)+abs(tagC6f(2:ncp,2:ncp,:)-jl);
0122 tagV(var==0)=jl;
0123 var=abs(tagC6f(2:ncp,1:nc,:)-jl)+abs(tagC6f(2:ncp,2:ncp,:)-jl+1);
0124 tagV(var==0)=-jl;
0125 end
0126
0127 %-------------------------------------------------------------------------------
0128 %- D.1) fourth step: make a list for each jl:
0129 % --> put back (reshape) tagU & tagV in 1 long vector
0130 % for each jl in 1:ydim
0131 % listU = find(abs(tagU)==jl); nU=length(listU); flag(1:nU)=sign(tagU)*1;
0132 % listV = find(abs(tagV)==jl); nV=length(listV); flag(1:nV)=sign(tagV)*2;
0133 tagU1d=reshape(tagU,[nc*nc*6 1]);
0134 tagV1d=reshape(tagV,[nc*nc*6 1]);
0135 for jl=1:ydim,
0136 [ijU]=find(abs(tagU1d)==jl);
0137 nU=length(ijU);
0138 [ijV]=find(abs(tagV1d)==jl);
0139 nV=length(ijV);
0140 %-
0141 n1=1; n2=nU;
0142 savFlag(n1:n2,jl)=sign(tagU1d(ijU));
0143 savIJuv(n1:n2,jl)=ijU;
0144 %-
0145 n1=nU+1; n2=nU+nV;
0146 savFlag(n1:n2,jl)=2*sign(tagV1d(ijV));
0147 savIJuv(n1:n2,jl)=ijV;
0148 savNpts(jl)=n2;
0149 end
0150 mx_pts=max(savNpts); [J]=find(savNpts==mx_pts);
0151 fprintf('Max Npts = %i (x %i) e.g. for jl= %i\n',mx_pts,length(J),J(1));
0152
0153 else
0154 %===============================================================================
0155 % Case with at least one Jump > 1 in adjacent tagC:
0156 % since one velocity pt can be listed in more than 1 bk-line, need to
0157 % process step (c) and (D) within jl loop.
0158 %-------------------------------------------------------------------------------
0159
0160 %tagU=zeros(nc,nc,6); tagV=tagU;
0161 for jl=1:ydim,
0162 tagU=zeros(nc,nc,6); tagV=tagU;
0163
0164 %- C.2) third step: Tag U & V point for this latitude yLat(jl):
0165 % tagU(i,j) = +jl if tagC6f(i,j) >= jl & tagC6f(i-1,j) < jl
0166 % tagU(i,j) = -jl if tagC6f(i,j) < jl & tagC6f(i-1,j) >= jl
0167 %
0168 % tagV(i,j) = +jl if tagC6f(i,j) >= jl & tagC6f(i,j-1) < jl
0169 % tagV(i,j) = -jl if tagC6f(i,j) < jl & tagC6f(i,j-1) >= jl
0170 %
0171 %- hard to include case == jl with this:
0172 %var=(tagC6f(1:nc,2:ncp,:)-jl).*(tagC6f(2:ncp,2:ncp,:)-jl);
0173 %tagU(var < 0)=jl;
0174 %tagU=tagU.*sign(tagC6f(2:ncp,2:ncp,:)-tagC6f(1:nc,2:ncp,:));
0175
0176 var=zeros(nc,nc,6);
0177 tmp=tagC6f(1:nc,2:ncp,:); var(tmp < jl)=1;
0178 tmp=tagC6f(2:ncp,2:ncp,:); var(tmp < jl)=0;
0179 tagU(var == 1)=jl;
0180
0181 var=zeros(nc,nc,6);
0182 tmp=tagC6f(2:ncp,2:ncp,:); var(tmp < jl)=1;
0183 tmp=tagC6f(1:nc,2:ncp,:); var(tmp < jl)=0;
0184 tagU(var == 1)=-jl;
0185
0186 var=zeros(nc,nc,6);
0187 tmp=tagC6f(2:ncp,1:nc,:); var(tmp < jl)=1;
0188 tmp=tagC6f(2:ncp,2:ncp,:); var(tmp < jl)=0;
0189 tagV(var == 1)=jl;
0190
0191 var=zeros(nc,nc,6);
0192 tmp=tagC6f(2:ncp,2:ncp,:); var(tmp < jl)=1;
0193 tmp=tagC6f(2:ncp,1:nc,:); var(tmp < jl)=0;
0194 tagV(var == 1)=-jl;
0195
0196 %-------------------------------------------------------------------------------
0197 %- D.1) fourth step: make a list for this latitude line yLat(jl):
0198 % --> put back (reshape) tagU & tagV in 1 long vector
0199 % for each jl in 1:ydim
0200 % listU = find(abs(tagU)==jl); nU=length(listU); flag(1:nU)=sign(tagU)*1;
0201 % listV = find(abs(tagV)==jl); nV=length(listV); flag(1:nV)=sign(tagV)*2;
0202 tagU1d=reshape(tagU,[nc*nc*6 1]);
0203 tagV1d=reshape(tagV,[nc*nc*6 1]);
0204 %---
0205 [ijU]=find(abs(tagU1d)==jl);
0206 nU=length(ijU);
0207 [ijV]=find(abs(tagV1d)==jl);
0208 nV=length(ijV);
0209 %-
0210 n1=1; n2=nU;
0211 savFlag(n1:n2,jl)=sign(tagU1d(ijU));
0212 savIJuv(n1:n2,jl)=ijU;
0213 %-
0214 n1=nU+1; n2=nU+nV;
0215 savFlag(n1:n2,jl)=2*sign(tagV1d(ijV));
0216 savIJuv(n1:n2,jl)=ijV;
0217 savNpts(jl)=n2;
0218
0219 end
0220
0221 end
0222
0223 %===============================================================================
0224 %- a quick plot to check:
0225 if kplot > 0,
0226 figure(1); clf; colormap('jet');
0227 ccB=[0 0]; shift=-1; cbV=1; AxBx=[-180 180 -90 90]; kEnv=0;
0228 var=tagC; titv='tagC';
0229 %var=tagU; titv='tagU';
0230 %var=tagV; titv='tagV';
0231 if strcmp(titv,'tagC'), ccB=[0 ydim];
0232 else
0233 if n2h == nc,
0234 var=reshape(permute(var,[1 3 2]),[n1h n2h]);
0235 else
0236 var=reshape(var,[n1h n2h]);
0237 end
0238 ccB=[-ydim ydim];
0239 end
0240 fprintf(' %s : min,Max= %9.5g , %9.5g \n',titv,min(var(:)),max(var(:)));
0241 if kplot == 1,
0242 grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
0243 else
0244 grph_CS_6t(var,ccB(1),ccB(2),nc,'tagC');
0245 end
0246 end
0247
0248 %--------------------------------------
0249 %- output : save bk-line into a matlab ".mat" file:
0250 %------------------------------
0251 if ydim > 1 & kwr == 1,
0252
0253 %- write smaller arrays : reduce size from nc*6 to Max(savNpts) :
0254 mxNpts=max(savNpts);
0255
0256 %- array to save in matlab "*.mat" file:
0257 % bkl_yLat = broke-line (bk-line) latidude-target vector
0258 % bkl_Npts = number of velocity points contributing to each bk-line
0259 % bkl_Flag = list to select sign and component (U or V) to use in bk-line
0260 % bkl_IJuv = index list of velocity points for each bk-line
0261 % bkl_tagC = 2-D map of grid-cell center tags, according to latitude band "yLat"
0262
0263 bkl_Flag=zeros(mxNpts,ydim);
0264 bkl_IJuv=zeros(mxNpts,ydim);
0265
0266 bkl_yLat=yLat ;
0267 bkl_Npts=savNpts ;
0268 bkl_Flag=savFlag(1:mxNpts,:);
0269 bkl_IJuv=savIJuv(1:mxNpts,:);
0270 if n2h == nc,
0271 bkl_tagC=permute(reshape(tagC,[nc 6 nc]),[1 3 2]);
0272 else
0273 bkl_tagC=reshape(tagC,[nc nc 6]);
0274 end
0275
0276 namf=['isoLat_',int2str(ydim),'_cs',int2str(nc)];
0277 save(namf,'bkl_yLat','bkl_Npts','bkl_Flag','bkl_IJuv','bkl_tagC');
0278 fprintf([' ==> Save %i lines on matlab file: ',namf,'\n'],ydim);
0279
0280 end
0281 %------------------------------
0282
0283 return
