Warning, /utils/matlab/cs_grid/bk_line/sep_API_basins.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
323aa12fe2 Jean*0001
0002 krd=1; kpr=1; kgr=0; kwr=1;
0003 %krd=0; kpr=0; kgr=1; kwr=0; % <- execute a 2nd time & draw some plot
0004
0005 if krd == 1,
0c7ba451ba Jean*0006 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0007 % or ncdf=0 to load MDS (binary) grid-files :
0008 ncdf=0;
0009 %rac='/home/jmc/grid_cs32/';
0010 rac='grid_files/';
0011
323aa12fe2 Jean*0012 %- load: bkl_Ylat, bkl_Npts, bkl_Flg, bkl_IJuv, bkl_Xsg, bkl_Ysg, bkl_Zon
0013 bk_lineF=[rac,'isoLat_cs32_59'];
0014 load(bk_lineF);
0015 ydim=length(bkl_Ylat); ydimC=ydim+1;
0016 fprintf([' load bk_line description from: ',bk_lineF,'.mat \n']);
0c7ba451ba Jean*0017
323aa12fe2 Jean*0018 %- load 'Nu','Nv','Iu','Iv':
0019 namf=[rac,'open_basins_section'];
0020 load(namf);
0c7ba451ba Jean*0021 fprintf([' and API separation from: ',namf,'.mat : O.K.\n']);
323aa12fe2 Jean*0022
0023 namf='maskC_bas.bin';
0024 fid=fopen([rac,namf],'r','b'); mskBasC=fread(fid,'real*4'); fclose(fid);
0025
f3103ba969 Jean*0026 G=load_grid(rac,0+ncdf);
0c7ba451ba Jean*0027 xcs=G.xC; ycs=G.yC; xcg=G.xG; ycg=G.yG; arc=G.rAc;
0028 ncx=G.dims(1); nc=G.dims(2); ncp=nc+1;
323aa12fe2 Jean*0029 %mskBasC=rdda([rac,namf],[ncx nc 3],1,'real*4','b');
0030 mskBasC=reshape(mskBasC,[ncx nc 3]);
0031
0032 xc1=reshape(xcs,1,ncx*nc);
0033 yc1=reshape(ycs,1,ncx*nc);
0034 %xg1=reshape(xcg,1,ncx*nc);
0035 %yg1=reshape(ycg,1,ncx*nc);
0036 %xg6=split_Z_cub(xcg);
0037 %yg6=split_Z_cub(ycg);
0038 xc6=split_C_cub(xcs,1);
0039 yc6=split_C_cub(ycs,1);
0040
0c7ba451ba Jean*0041 hFacC=G.hFacC; hFacW=G.hFacW; hFacS=G.hFacS;
323aa12fe2 Jean*0042 mskC=min(1,hFacC(:,:,1)); mskC=ceil(mskC);
0043 mskW=min(1,hFacW(:,:,1)); mskW=ceil(mskW);
0044 mskS=min(1,hFacS(:,:,1)); mskS=ceil(mskS);
0045 mskW=reshape(mskW,1,ncx*nc);
0046 mskS=reshape(mskS,1,ncx*nc);
0047
0048 end
0049
f3103ba969 Jean*0050 if kpr ==1,
323aa12fe2 Jean*0051 [xdum,xPA,yPA,xAI,yAI,xIP,yIP]=line_sep(0);
0c7ba451ba Jean*0052
0053 %- define Bas-Sep point (as U,V point) for each Lat-band
0054 % the list of U,V point
323aa12fe2 Jean*0055 MxSiz=ceil(nc*nc/2/ydimC);
0056 nbpSep=zeros(3,ydimC);
0057 ijSep=zeros(3,ydimC,MxSiz); typSep=zeros(3,ydimC,MxSiz);
0c7ba451ba Jean*0058 n2s=zeros(3,1); jl2s=zeros(3,MxSiz);
323aa12fe2 Jean*0059 xc2s=zeros(3,MxSiz,2); yc2s=zeros(3,MxSiz,2);
0060 for b=1:3,
0061 fprintf(' b= %i ; Nu,Nv= %3i %3i \n',b,Nu(b),Nv(b));
f3103ba969 Jean*0062 %-
0063 nU=Nu(b); I=Iu(b,1:nU); xSp=zeros(2,nU);
323aa12fe2 Jean*0064 i0=rem(I-1,ncx);j1=1+fix((I-1)/ncx);
0065 i1=1+rem(i0,nc); k1=1+fix(i0/nc);
0066 for i=1:nU,
0067 jl=1+bkl_Zon(I(i));
0068 nn=1+nbpSep(b,jl);
0069 nbpSep(b,jl)=nn;
0070 ijSep(b,jl,nn)=I(i);
0071 typSep(b,jl,nn)=1;
0072 if i1(i) == 1,
0073 xx=xc6(i1(i),1+j1(i),k1(i)); yy=yc6(i1(i),1+j1(i),k1(i));
0074 [Im]=find( xc1 == xx & yc1 == yy );
0075 if length(Im) ~= 1, fprintf(' stop A \n'); return ; end
0076 else Im=I(i)-1; end
0077 jj=1+bkl_Zon(Im);
0078 if jj ~= jl,
0079 nn=1+nbpSep(b,jj);
0080 nbpSep(b,jj)=nn;
0081 ijSep(b,jj,nn)=I(i);
0082 typSep(b,jj,nn)=-1;
0083 if (mskW(I(i)) == 1),
0084 nn=n2s(b)+1; n2s(b)=nn;
0085 jl2s(b,nn)=max(jj,jl);
0086 xc2s(b,nn,1)=xc1(Im); yc2s(b,nn,1)=yc1(Im);
0087 xc2s(b,nn,2)=xc1(I(i)); yc2s(b,nn,2)=yc1(I(i));
0088 end
0089 end
f3103ba969 Jean*0090 xSp(1,i)=xc6(1+i1(i),1+j1(i),k1(i));
0091 xSp(2,i)=xc6( i1(i), 1+j1(i),k1(i));
323aa12fe2 Jean*0092 ySepU(b,i)=(yc6(1+i1(i),1+j1(i),k1(i))+yc6(i1(i),1+j1(i),k1(i)))/2;
0093 end
f3103ba969 Jean*0094 var=reshape(xcg,[1 ncx*nc]);
0095 xSp(1,:)=var(I)-180+rem(xSp(1,:)-var(I)+3*180,360);
0096 xSp(2,:)=var(I)-180+rem(xSp(2,:)-var(I)+3*180,360);
0097 xSepU(b,1:nU)=(xSp(1,:)+xSp(2,:))/2;
0098 %-
0099 nV=Nv(b); J=Iv(b,1:nV); xSp=zeros(2,nV);
323aa12fe2 Jean*0100 i0=rem(J-1,ncx);j1=1+fix((J-1)/ncx);
0101 i1=1+rem(i0,nc); k1=1+fix(i0/nc);
0102 for i=1:nV,
0103 jl=1+bkl_Zon(J(i));
0104 nn=1+nbpSep(b,jl);
0105 nbpSep(b,jl)=nn;
0106 ijSep(b,jl,nn)=J(i);
0107 typSep(b,jl,nn)=2;
0108 if j1(i) == 1,
0109 xx=xc6(1+i1(i),j1(i),k1(i)); yy=yc6(1+i1(i),j1(i),k1(i));
0110 [Jm]=find( xc1 == xx & yc1 == yy );
0111 if length(Jm) ~= 1, fprintf(' stop B \n'); return ; end
0112 else Jm=J(i)-ncx; end
0113 jj=1+bkl_Zon(Jm);
0114 if jj ~= jl,
0115 nn=1+nbpSep(b,jj);
0116 nbpSep(b,jj)=nn;
0117 ijSep(b,jj,nn)=J(i);
0118 typSep(b,jj,nn)=-2;
0119 if (mskS(J(i)) == 1),
0120 nn=n2s(b)+1; n2s(b)=nn;
0121 jl2s(b,nn)=max(jj,jl);
0122 xc2s(b,nn,1)=xc1(Jm); yc2s(b,nn,1)=yc1(Jm);
0123 xc2s(b,nn,2)=xc1(J(i)); yc2s(b,nn,2)=yc1(J(i));
0124 end
0125 end
f3103ba969 Jean*0126 xSp(1,i)=xc6(1+i1(i),1+j1(i),k1(i));
0127 xSp(2,i)=xc6(1+i1(i), j1(i), k1(i));
323aa12fe2 Jean*0128 ySepV(b,i)=(yc6(1+i1(i),1+j1(i),k1(i))+yc6(1+i1(i),j1(i),k1(i)))/2;
0129 end
f3103ba969 Jean*0130 xSp(1,:)=var(J)-180+rem(xSp(1,:)-var(J)+3*180,360);
0131 xSp(2,:)=var(J)-180+rem(xSp(2,:)-var(J)+3*180,360);
0132 xSepV(b,1:nV)=(xSp(1,:)+xSp(2,:))/2;
323aa12fe2 Jean*0133 end
0134 NbpSep=max(nbpSep');
0135 fprintf(' NbpSep = %2i %2i %2i \n',NbpSep);
0136 fprintf(' Dbl Sep: %2i %2i %2i \n',n2s);
0137
0138 end
0139
f3103ba969 Jean*0140 if kwr ==1,
323aa12fe2 Jean*0141 %- reduce size :
0142 mxSiz=max(NbpSep);
0143 ij_Sep=ijSep(:,:,1:mxSiz); tp_Sep=typSep(:,:,1:mxSiz); np_Sep=nbpSep;
0144
0145 namf=['sepBas_cs',int2str(nc),'_',int2str(ydimC)];
0146 save(namf,'np_Sep','ij_Sep','tp_Sep');
0147 fprintf(['write Basin Separation on file:',namf,'.mat : O.K. \n']);
0148 end
0149
0c7ba451ba Jean*0150 if kgr ==1,
0151 figure(1); clf;
f3103ba969 Jean*0152 shift=-1; cbV=1; ccB=[0 0]; AxBx=[-182 182 -90 90];
0153 for bp=1:3, subplot(310+bp);
0154 %for bp=2:2,
0155 %AxBx=[80 181 30 91];
0156 %AxBx=[80 181 -91 -30];
0c7ba451ba Jean*0157 var=mskBasC(:,:,1)+2*mskBasC(:,:,2)+3*mskBasC(:,:,3);
323aa12fe2 Jean*0158 ccB=[-1 5]; var(find(mskC==0))=NaN;
0159 grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
0160 hold on ;
f3103ba969 Jean*0161 [L0]=plot(xPA,yPA,'*b-',xAI,yAI,'*b-',xIP,yIP,'*b-'); set(L0,'MarkerSize',6);
323aa12fe2 Jean*0162 for b=bp:bp,
0163 xloc=squeeze(xc2s(b,:,:)); yloc=squeeze(yc2s(b,:,:));
0164 for i=1:n2s(b), jl=jl2s(b,i)-1;
0165 ie=bkl_Npts(jl)+1;
0166 [L1]=plot(bkl_Xsg(1:ie,jl),bkl_Ysg(1:ie,jl),'r-');
0167 %set(L1,'Color',[.9 .2 .0],'LineWidth',1); % set(L1,'LineStyle','-');
f3103ba969 Jean*0168 [L2]=plot(xloc(i,:),yloc(i,:),'xw-');
323aa12fe2 Jean*0169 set(L2,'LineWidth',2);
0170 end
f3103ba969 Jean*0171 [L3]=plot(xSepU(b,1:Nu(b)),ySepU(b,1:Nu(b)),'m*'); set(L3,'MarkerSize',6);
0172 [L4]=plot(xSepV(b,1:Nv(b)),ySepV(b,1:Nv(b)),'mo'); set(L4,'MarkerSize',6);
323aa12fe2 Jean*0173 end
0174 hold off
0175 end
0176 end
0177
0178 return
