cs_grid/bk_line/gen_bk_Zon.m

323aa12fe2 Jean*0001 krd=1; kpr=1; kgr=0; kwr=1;
                0002 % krd=0; kpr=0; kgr=1; kwr=0; % <- execute a 2nd time & draw some plot
                0003 ijprt=65;
                0004
                0005 if krd == 1,
                0006
0c7ba451ba Jean*0007 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
                0008 %   or ncdf=0 to load MDS (binary) grid-files :
                0009  ncdf=0;
                0010
323aa12fe2 Jean*0011 %rac='/home/jmc/grid_cs32/';
0c7ba451ba Jean*0012  rac='grid_files/';
323aa12fe2 Jean*0013  bk_lineF=[rac,'isoLat_cube32_59'];
                0014
                0015 %- load broken lines :
                0016 % bkl_Ylat,bkl_Npts,bkl_Flg,bkl_Iuv,bkl_Juv,bkl_Xsg,bkl_Ysg
                0017  load(bk_lineF);
                0018
0c7ba451ba Jean*0019 %- load grid-files :
                0020  G=load_grid(rac,10+ncdf);
                0021  xcs=G.xC; ycs=G.yC; xcg=G.xG; ycg=G.yG; arc=G.rAc;
                0022  ncx=G.dims(1); nc=G.dims(2);
                0023
                0024 %load_cs;
                0025 %ncx=size(xcs,1); nc=size(xcs,2);
323aa12fe2 Jean*0026
                0027 %xc6=split_C_cub(xcs,1);
                0028 %yc6=split_C_cub(ycs,1);
                0029 %xg6=split_Z_cub(xcg);
                0030 %yg6=split_Z_cub(ycg);
                0031  xc1=reshape(xcs,ncx*nc,1);
                0032  yc1=reshape(ycs,ncx*nc,1);
                0033
                0034 %rac='res_r1/';
                0035 %hFacC=rdmds([rac,'hFacC']);
                0036 %mskC=min(1,hFacC); mskC=ceil(mskC);
                0037
                0038  fprintf([' load bk_line description from file=',bk_lineF,'\n']);
                0039 end
                0040
                0041 ydim=length(bkl_Ylat); ylat=zeros(ydim+2,1);
                0042 ylat(2:ydim+1)=bkl_Ylat; ylat(1)=-90; ylat(ydim+2)=90;
                0043 bkl_IJuv=bkl_Iuv+ncx*(bkl_Juv-1);
                0044 Yg1=zeros(ydim,1); Yg2=Yg1;
                0045 for jl=1:ydim,
                0046   ie=bkl_Npts(jl)+1;
                0047   Yg1(jl)=min(bkl_Ysg(1:ie,jl)); Yg2(jl)=max(bkl_Ysg(1:ie,jl));
                0048 end
                0049
                0050 if kpr == 1,
                0051
                0052  bkl_Zon=-ones(ncx*nc,1);
                0053  for ij=1:ncx*nc,
                0054   xx=xc1(ij); yy=yc1(ij);
0c7ba451ba Jean*0055   [J2]=find(Yg1 >= yy);
323aa12fe2 Jean*0056   if length(J2)==0, jl2=ydim+1; else jl2=min(J2); end
0c7ba451ba Jean*0057   [J1]=find(Yg2 <= yy);
323aa12fe2 Jean*0058   if length(J1)==0, jl1=0; else jl1=max(J1); end
                0059   fprintf('ij,xx,yy= %4i %6.1f %6.1f ; jl1,jl2= %2i %2i',ij, xx, yy, jl1,jl2);
                0060   %- solve for intersection bkl_X,Y with x=xx :
                0061   jL1=-1; jL2=-1; yXX=zeros(1,1+jl2-jl1);
                0062   for jl=jl1:jl2
                0063    if ij==ijprt, fprintf('\n  jl= %i ',jl);end
0c7ba451ba Jean*0064    l=1+jl-jl1;
323aa12fe2 Jean*0065    if jl == 0, yXX(l)=-90 ;
0c7ba451ba Jean*0066    elseif jl == ydim+1, yXX(l)=90;
323aa12fe2 Jean*0067    else
                0068     ie=bkl_Npts(jl)+1;
                0069     xloc=bkl_Xsg(1:ie,jl); yloc=bkl_Ysg(1:ie,jl); yloc(ie+1)=yloc(1);
0c7ba451ba Jean*0070     xloc2=xloc; xloc2(1:ie-1)=xloc(2:ie); xloc2(ie)=xloc(1)+360;
323aa12fe2 Jean*0071     if ij==ijprt, fprintf('; xloc(1)= %6.1f ',xloc(1)); end
                0072     [I]=find( xloc <= xx & xx < xloc2 );
                0073     if length(I) ~= 1, fprintf(' Stop A \n');return; end
                0074     if ij==ijprt, fprintf('; x,y_loc= %6.1f %6.1f %5.1f %5.1f ', ...
                0075                             xloc(I),xloc2(I),yloc(I),yloc(I+1)); end
                0076     if xloc2(I) == xloc(I),
                0077      yXX(l)=yloc(I);
                0078     else
                0079      yXX(l)=yloc(I)+(yloc(I+1)-yloc(I))*(xx-xloc(I))/(xloc2(I)-xloc(I));
                0080     end
                0081    end
                0082    if ij==ijprt, fprintf('; yXX= %6.1f',yXX(l)); end
                0083    if l>1 & yXX(l) < yXX(l-1), fprintf(' Stop B \n');return; end
                0084    if yXX(l) <= yy, jL1=jl; end
                0085    if yXX(l) > yy & jL2==-1, jL2=jl; end
                0086   end
0c7ba451ba Jean*0087   if jL2 == jL1+1, bkl_Zon(ij)=jL1;
323aa12fe2 Jean*0088    fprintf('; jL1= %2i\n',jL1);
                0089   else fprintf(' Stop C : %i %i %6.1f %6.1f\n',jL1,jL2,xx,yy);
                0090    fprintf(' yXX :');fprintf(' %9.5f',yXX(1:1+jl2-jl1)); fprintf('\n');
0c7ba451ba Jean*0091    return;
323aa12fe2 Jean*0092   end
                0093  end % <- ij loop
                0094
                0095 end
                0096
                0097 if kwr==1,
                0098  namf=['isoLat_cs',int2str(nc),'_',int2str(ydim)];
                0099 %- save to matlab file:
                0100 % bkl_Ylat, bkl_Npts, bkl_Flg, bkl_IJuv, bkl_Xsg, bkl_Ysg, bkl_Zon
                0101  save(namf,'bkl_Ylat','bkl_Npts','bkl_Flg','bkl_IJuv','bkl_Xsg','bkl_Ysg','bkl_Zon');
                0102  fprintf([' ==> Save %i lines on matlab file: ',namf,'\n'],ydim);
                0103 end
                0104
                0105 if kgr==1,
                0106  figure(1); clf;
040160fb38 Andr*0107  shift=-1; cbV=2; ccB=[0 0]; AxBx=[-180 180 -90 90];
323aa12fe2 Jean*0108 %var=zeros(ncx,nc) ; var(1+rem(ijprt-1,ncx),1+fix((ijprt-1)/ncx))=1 ; ccB=[-1 1];
                0109  var=reshape(bkl_Zon,ncx,nc) ; ccB=[-1 50];
                0110  grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
                0111  hold on ;
                0112  for jl=jl1:jl2,
                0113   ie=bkl_Npts(jl)+1;
                0114   [L]=plot(bkl_Xsg(1:ie,jl),bkl_Ysg(1:ie,jl),'r-');
                0115   set(L,'Color',[0 0 0],'LineWidth',2); % set(L,'LineStyle','-');
                0116  end
                0117  plot(xx*ones(1,1+jl2-jl1),yXX,'ro-');
                0118  hold off
                0119 end