cs_grid/latloncap/add_Angle2grid.m

36056b62bf Jean*0001
                0002 %-- Add angle (Cos & Sin) to grid files
                0003
                0004 %nr = 360; ng = 120; nb = 120; tnx=120; tny=120; %- polar-cap grid
                0005 rgbDim=[360 120 120];
                0006 %- set all 6 faces dimensions
                0007  nr=rgbDim(1); ng=rgbDim(2); nb=rgbDim(3);
                0008  nf=ones(6,2);
                0009  nf(1,:)=[nb nr]; nf(2,:)=[ng nr]; nf(3,:)=[ng nb];
                0010  nf(4,:)=[nr nb]; nf(5,:)=[nr ng]; nf(6,:)=[nb ng];
                0011  fdim=prod(nf,2); fd2= cumsum(fdim); fd1=fd2-fdim+1;
                0012
                0013 %- compact format:
                0014 nx=gcd(gcd(rgbDim(1),rgbDim(2)),rgbDim(3));
                0015 ny=sum(fdim)/nx;
                0016
                0017 %--------
                0018 %rDir='../input/';
                0019 rDir='./';
                0020 inpName='face.mitgrid';
                0021 kwr=1; outpName='llc_120x360';
                0022
                0023 for n=1:6,
                0024 %-- read :
                0025  p=strfind(inpName,'.');
                0026  if isempty(p),
                0027    namF=sprintf([rDir,inpName,'.face%3.3i.bin'],n);
                0028  else
                0029    namF=sprintf([rDir,inpName(1:4),'%3.3i',inpName(5:end)],n);
                0030  end
                0031  fid=fopen(namF,'r','b');
                0032  vv1=fread(fid,'real*8');
                0033  fclose(fid);
                0034  s=size(vv1,1); k1=s/(nf(n,1)+1)/(nf(n,2)+1);
                0035  fprintf(['read: ',namF,' : size: %i (%ix%ix%i)'],s,nf(n,1)+1,nf(n,2)+1,k1);
                0036  vv1=reshape(vv1,[nf(n,1)+1 nf(n,2)+1 k1]);
                0037 %- XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG:
                0038 %------------------------------------------------------------------
                0039 %-- save each face:
                0040  nvar=['vvf',int2str(n),'=vv1;']; fprintf(' ; %s \n',nvar);
                0041  eval(nvar);
                0042 end
                0043 titv={'XC','YC','DXF','DYF','RA','XG','YG','DXV','DYU','RAZ','DXC','DYC', ...
                0044       'RAW','RAS','DXG','DYG'};
                0045
                0046 jlc=330; ilc=1+nr-jlc;
                0047
                0048 %- XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG:
                0049 %   1  2  3   4   5  6  7  8   9  10  11  12  13  14  15  16
                0050 gpos=[3 3 3   3   3  0  0  0   0   0   1   2   1   2   2   1];
                0051 % gpos: 0 = grid-cell corner ; 1 = U point ; 2 = V point ; 3 = Tracer point
                0052
                0053 %----------------------------------------------------------------------------
                0054
                0055 %- compute the analytic grid :
                0056 xc1=vvf1(1:120,1:jlc,1); xc=mean(xc1,2);
                0057 yc1=vvf1(1:120,1:jlc,2); yc=mean(yc1,1);
                0058 xg1=vvf1(1:121,1:jlc+1,6); xg=mean(xg1,2);
                0059 yg1=vvf1(1:121,1:jlc+1,7); yg=mean(yg1,1);
                0060
                0061 %- check symetry:
                0062 jeq=find(abs(yg) < 1.e-4);
                0063 dxEq=vvf1(1:120,jeq,15); Radius=2*sum(dxEq)/pi;
                0064 %fprintf(' Radius= %f ',Radius);
                0065 %Radius=round(Radius)*1; dxDeg=0.75;
                0066 fprintf(', Radius= %f\n',Radius);
                0067 rad=pi/180;
                0068
                0069 k2=k1+2;
                0070 for n=1:6,
                0071  nvar=['vv1=vvf',int2str(n),';']; %fprintf(' ; %s \n',nvar);
                0072  eval(nvar);
                0073  i2=nf(n,1); j2=nf(n,2); ip=i2+1; jp=j2+1;
                0074  vv2=zeros(ip,jp,k2); vv2(:,:,[1:k1])=vv1;
                0075  if n < 6,
                0076   yg=vv1(:,:,7);
                0077   dxg=vv1(1:i2,1:jp,15);
                0078   dyg=vv1(1:ip,1:j2,16);
                0079 % Build purely zonal flow from a stream function, psi = r_earth * latitude.
                0080   psi=-Radius*yg*rad;
                0081   uZ=psi([1:ip],[1:j2],:)-psi([1:ip],[2:jp],:);
                0082   vZ=psi([2:ip],[1:jp],:)-psi([1:i2],[1:jp],:);
                0083   uZ=uZ./dyg;
                0084   vZ=vZ./dxg;
                0085 % Put constructed zonal wind at cell center.
                0086   uZc=(uZ(1:i2,:,:)+uZ(2:ip,:,:))/2;
                0087   vZc=(vZ(:,1:j2,:)+vZ(:,2:jp,:))/2;
                0088 % Calculate angles.
                0089   norm=sqrt(uZc.*uZc+vZc.*vZc);
                0090   AngleCS =  uZc./norm;
                0091   AngleSN = -vZc./norm;
                0092   vv2(1:i2,1:j2,k1+1)=AngleCS;
                0093   vv2(1:i2,1:j2,k1+2)=AngleSN;
                0094  else
                0095   AngleCS =  ones(i2,j2);
                0096   AngleSN =  zeros(i2,j2);
                0097  end
                0098 %-- save to Structure:
                0099  nvar=sprintf('csF.f%3.3i=AngleCS;',n); eval(nvar)
                0100  nvar=sprintf('snF.f%3.3i=AngleSN;',n); eval(nvar)
                0101
                0102 %- write to file:
                0103  if kwr == 1,
                0104 %-- write :
                0105   namF=sprintf([rDir,outpName,'.face%3.3i.bin'],n);
                0106   fprintf([': write to ',namF,' : size: (%ix%ix%i)\n'],nf(n,1)+1,nf(n,2)+1,k2);
                0107   fid=fopen(namF,'w','b');
                0108   fwrite(fid,vv2,'real*8');
                0109   fclose(fid);
                0110  end
                0111
                0112 end
                0113
                0114 nfg=1; k=0; ccB=[-1 1]*1.5;
                0115 plot_faces(nfg,csF,k,ccB,rgbDim);
                0116 nfg=2;
                0117 plot_faces(nfg,snF,k,ccB,rgbDim);
                0118 %----------------------------------------------------------------------------
                0119 return