cs_grid/latloncap/read_llc_fkij.m

ee6c6a54c1 Dimi*0001 function fld=read_llc_fkij(fnam,nx,face,kx,ix,jx,prec)
                0002
4e55f37103 Dimi*0003 % Function fld=read_llc_fkij(fnam,nx,face,kx,ix,jx,prec)
ee6c6a54c1 Dimi*0004 % read in specific face and indices for llc configuration
                0005 %
                0006 % INPUTS
                0007 % fnam  input path and file name
                0008 % nx    tile dimension (default 270)
                0009 % face  face number, 1 to 5, to read (default 1)
                0010 % kx    vertical indices to read, e.g., 1:50 (default 1)
                0011 % ix    i-indices, 1 to nx, to read (default 1:nx)
                0012 % jx    j-indices, 1 to 3*nx, to read (default 1:nx or 1:3*nx)
                0013 % prec  numeric precision (see fread; default 'real*4')
                0014 %
                0015 % OUTPUTS
                0016 % fld   output array of dimension length(ix)*length(jx)*length(kx)
                0017 %
e6b73c8b78 Dimi*0018 % NOTES
72bbcd60d4 Dimi*0019 %  For faces 4 and 5, model fields defined at the edges are:
                0020 %   u_East (i,j) =   v_read_llc(i,j)
                0021 %   v_North(i,j) = - u_read_llc(i,j-1)   <= note “j-1"
                0022 %   dyG(i,j)     = dxG_read_llc(i,j)
                0023 %   dxG(i,j)     = dyG_read_llc(i,j-1)   <= again note “j-1"
e6b73c8b78 Dimi*0024 %
ee6c6a54c1 Dimi*0025 % SEE ALSO
                0026 % quilread_llc
                0027
                0028 if nargin < 1, error('please specify input file name'); end
                0029 if nargin < 2, nx=270; end
                0030 if nargin < 3, face=1; end
                0031 if nargin < 4, kx=1; end
                0032 if nargin < 5, ix=1:nx; end
                0033 if nargin < 6,
                0034     jx=1:(3*nx);
                0035     if face==3, jx=1:nx; end
                0036 end
                0037 if nargin < 7, prec='real*4'; end
                0038
                0039 % reverse jx if reading faces 4 and 5
                0040 if face > 3
                0041     jx=3*nx+1-jx(end:-1:1);
                0042 end
                0043
                0044 fid=fopen(fnam,'r','ieee-be');
                0045
                0046 switch prec
f91ed0bfc5 Dimi*0047  case {'integer*1'}
                0048   prec='int8';
                0049  case {'integer*2'}
                0050   prec='int16';
                0051  case {'integer*4'}
                0052   prec='int32';
                0053  case {'real*4','float32'}
                0054   prec='single';
                0055  case {'integer*8'}
                0056   prec='int64';
                0057  case {'real*8','float64'}
                0058   prec='double';
                0059 end
                0060
                0061 switch prec
                0062  case {'int8'}
ee6c6a54c1 Dimi*0063   preclength=1;
f91ed0bfc5 Dimi*0064  case {'int16','uint16'}
ee6c6a54c1 Dimi*0065   preclength=2;
f91ed0bfc5 Dimi*0066  case {'int32','uint32','single'}
ee6c6a54c1 Dimi*0067   preclength=4;
f91ed0bfc5 Dimi*0068  case {'int64','uint64','double'}
ee6c6a54c1 Dimi*0069   preclength=8;
                0070 end
                0071
f91ed0bfc5 Dimi*0072 fld=zeros(length(ix),length(jx),length(kx),prec);
ee6c6a54c1 Dimi*0073
                0074 for k=1:length(kx)
                0075     skip=(kx(k)-1)*nx*nx*13; % numbers to skip to reach vertical level kx(k)
                0076     switch face
                0077       case {1,2,3}
                0078         skip=skip + ...
                0079              (face-1)*3*nx*nx + ...     % add numbers to reach face
                0080              (jx(1)-1)*nx + ...         % add numbers to reach row jx(1)
                0081              ix(1)-1;                   % add numbers to reach column ix(1)
                0082         if(fseek(fid,skip*preclength,'bof')<0)
                0083             error(['past end of file']);
                0084         end
                0085         tmp=fread(fid,[length(ix),length(jx)], ...
                0086                   [int2str(length(ix)) '*' prec], ... % repetition factor
                0087                   (nx-length(ix))*preclength);        % skip argument
                0088       case {4,5}
                0089         oskip=skip;
                0090         skip=skip + ...
                0091              ((face-1)*3-2)*nx*nx + ... % add numbers to reach face
                0092              (ix(1)-1)*3*nx + ...       % add numbers to reach row ix(1)
                0093              jx(1)-1;                   % add numbers to reach column jx(1)
                0094         if(fseek(fid,skip*preclength,'bof')<0)
                0095             error(['past end of file']);
                0096         end
                0097         tmp=rot90(fread(fid,[length(jx),length(ix)], ...
                0098                         [int2str(length(jx)) '*' prec], ... % repetition factor
                0099                         (3*nx-length(jx))*preclength), ...  % skip argument
                0100                   3);
                0101     end
                0102     fld(:,:,k)=tmp;
                0103 end
                0104
                0105 fid=fclose(fid);