detection/read_wrfout_sel.m

   1 function h=read_wrfout_sel(files,vars,nsamples)
   2 % h=read_wrfout_sel(files,vars,nsamples)
   3 % read selected time levels from wrfouts
   4 % in
   5 %     files     cell array of file names
   6 %     files     cell array of variable names; if followed by 'sparse' will be read as such
   7 %     nsamples  number of timestep samples to read, if missing read all
   8 % out
   9 %     h         structure with the selected variables
  10
  11 for i=1:length(files)
  12     [s,dims]=nc2struct(files{i},{'Times'},{});
  13     n(i)=dims.times(2);
  14 end
  15 nvar=0;
  16 lastvar=0;
  17 for i=1:length(vars)
  18     switch vars{i}
  19     case 'sparse'
  20         if ~lastvar,
  21             error('sparse flag must follow a variable name');
  22         end
  23         varsparse(nvar)=1;
  24         lastvar=0;
  25     otherwise
  26         lastvar=1;
  27         nvar=nvar+1;
  28         varname{nvar}=vars{i};
  29         spfield{nvar}=[lower(varname{nvar}),'_sparse'];
  30         varsparse(nvar)=0;
  31     end
  32 end
  33
  34 ntot=sum(n);
  35 if ~exist('nsamples','var') | nsamples == 0,
  36     step=1;
  37 else
  38     step = floor((ntot+nsamples-1)/nsamples)
  39 end
  40 start=[1,cumsum(n)+1]
  41 idx=mod(ntot-1,step)+1:step:ntot;
  42 nidx=length(idx);
  43 for jx=1:nidx
  44     j=idx(jx);
  45     k=find(j>=start);k=k(end);
  46     loc=j-start(k)+1;
  47     fprintf('reading step %i as %i from file %i into %i\n',j,loc,k,jx)
  48     if loc <=0 | k > length(files), error('bad index'), end
  49     [f,dims]=nc2struct(files{k},varname,{},loc);
  50     if ~exist('h','var'),
  51         for jj=1:length(varname)
  52             field=lower(varname{jj});
  53             d=dims.(field);
  54             if varsparse(jj),
  55                 h.(spfield{jj})=cell(1,nidx);
  56             else
  57                 h.(field)=zeros([d(1:end-1),nidx]);
  58             end
  59         end
  60     end
  61     for jj=1:length(varname)
  62         if varsparse(jj),
  63             h.(spfield{jj}){jx}=sparse(f.(field));
  64         else
  65             field=lower(varname{jj});
  66             n=length(dims.(field));
  67             switch n
  68             case 2
  69                 h.(field)(:,jx)=f.(field);
  70             case 3
  71                 h.(field)(:,:,jx)=f.(field);
  72             case 4
  73                 h.(field)(:,:,:,jx)=f.(field);
  74             otherwise
  75                 error(['unsupported number of dimensions ',num2str(n)]);
  76             end
  77         end
  78     end
  79 end
  80
  81