/[MITgcm]/MITgcm/utils/matlab/rdmnc.m
ViewVC logotype

Diff of /MITgcm/utils/matlab/rdmnc.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.18 by mlosch, Fri Nov 13 07:57:16 2009 UTC revision 1.19 by mlosch, Tue Nov 24 13:09:16 2009 UTC
# Line 148  function [i0,j0,fn] = findTileOffset(S); Line 148  function [i0,j0,fn] = findTileOffset(S);
148    
149  function [S] = rdmnc_local(nc,varlist,iters,S,dBug)  function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
150    
151      fiter = nc{'iter'}(:);                               % File iterations present    fiter = nc{'iter'}(:);                               % File iterations present
152      if isempty(fiter), fiter = nc{'T'}(:); end    if isempty(fiter), fiter = nc{'T'}(:); end
153      if isinf(iters); iters = fiter(end); end    if isinf(iters); iters = fiter(end); end
154      if isnan(iters); iters = fiter; end    if isnan(iters); iters = fiter; end
155      [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index    [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
156      dii = dii(find(dii ~= 0));                           % Data interation index    dii = dii(find(dii ~= 0));                           % Data interation index
157      if dBug > 0,    if dBug > 0,
158        fprintf(' ; fii='); fprintf(' %i',fii);      fprintf(' ; fii='); fprintf(' %i',fii);
159        fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');      fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
160      end    end
161            
162      % No variables specified? Default to all    % No variables specified? Default to all
163      if isempty(varlist), varlist=ncnames(var(nc)); end    if isempty(varlist), varlist=ncnames(var(nc)); end
164          
165      % Attributes for structure    % Attributes for structure
166      if iters>0; S.iters_from_file=iters; end    if iters>0; S.iters_from_file=iters; end
167      S.attributes.global=read_att(nc);    S.attributes.global=read_att(nc);
168      [pstr,netcdf_fname,ext] = fileparts(name(nc));    [pstr,netcdf_fname,ext] = fileparts(name(nc));
169      if strcmp(netcdf_fname(end-3:end),'glob')    if strcmp(netcdf_fname(end-3:end),'glob')
170        % assume it is a global file produced by gluemnc and change some      % assume it is a global file produced by gluemnc and change some
171        % attributes      % attributes
172        S.attributes.global.sNx = S.attributes.global.Nx;      S.attributes.global.sNx = S.attributes.global.Nx;
173        S.attributes.global.sNy = S.attributes.global.Ny;      S.attributes.global.sNy = S.attributes.global.Ny;
174        S.attributes.global.nPx = 1;      S.attributes.global.nPx = 1;
175        S.attributes.global.nSx = 1;      S.attributes.global.nSx = 1;
176        S.attributes.global.nPy = 1;      S.attributes.global.nPy = 1;
177        S.attributes.global.nSy = 1;      S.attributes.global.nSy = 1;
178        S.attributes.global.tile_number = 1;      S.attributes.global.tile_number = 1;
179        S.attributes.global.nco_openmp_thread_number = 1;      S.attributes.global.nco_openmp_thread_number = 1;
180      end    end
181            
182          % Read variable data    % Read variable data
183          for ivar=1:size(varlist,2)    for ivar=1:size(varlist,2)
184                
185          cvar=char(varlist{ivar});      cvar=char(varlist{ivar});
186          if isempty(nc{cvar})      if isempty(nc{cvar})
187              disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);        disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
188              continue        continue
189          end      end
190                
191          dims = ncnames(dim(nc{cvar}));        % Dimensions      dims = ncnames(dim(nc{cvar}));        % Dimensions
192          sizVar = size(nc{cvar}); nDims=length(sizVar);      sizVar = size(nc{cvar}); nDims=length(sizVar);
193          if dims{1} == 'T'      if dims{1} == 'T'
194              if isempty(find(fii)), error('Iters not found'); end        if isempty(find(fii)), error('Iters not found'); end
195              it = length(dims);        it = length(dims);
196              tmpdata = nc{cvar}(fii,:);        tmpdata = nc{cvar}(fii,:);
197  %-      leading unity dimensions get lost; add them back:        % leading unity dimensions get lost; add them back:
198              tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);        tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
199          else      else
200              it = 0;        it = 0;
201              tmpdata = nc{cvar}(:);        tmpdata = nc{cvar}(:);
202  %-      leading unity dimensions get lost; add them back:        % leading unity dimensions get lost; add them back:
203              tmpdata=reshape(tmpdata,sizVar);        tmpdata=reshape(tmpdata,sizVar);
204          end      end
205                
206          if dBug > 1,      if dBug > 1,
207            fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));        fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
208            fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);        fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
209          end      end
210          if length(dims) > 1,      if length(dims) > 1,
211            tmpdata=permute(tmpdata,[nDims:-1:1]);        tmpdata=permute(tmpdata,[nDims:-1:1]);
212          end      end
213          if dBug > 1,      if dBug > 1,
214  %         fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');  %         fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
215          end      end
216          [ni nj nk nm nn no np]=size(tmpdata);      [ni nj nk nm nn no np]=size(tmpdata);
217                  
218          [i0,j0,fn]=findTileOffset(S);      [i0,j0,fn]=findTileOffset(S);
219          cdim=dims{end}; if cdim(1)~='X'; i0=0; end      cdim=dims{end}; if cdim(1)~='X'; i0=0; end
220          cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end      cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
221          if length(dims)>1;      if length(dims)>1;
222              cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end        cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
223          else      else
224              j0=0;        j0=0;
225          end      end
226          if dBug > 1,      if dBug > 1,
227            fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);        fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
228          end      end
229                
230          Sstr = '';      Sstr = '';
231          for istr = 1:max(nDims,length(dims)),      for istr = 1:max(nDims,length(dims)),
232              if     istr == it,  Sstr = [Sstr,'dii,'];        if     istr == it,  Sstr = [Sstr,'dii,'];
233              elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];        elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
234              elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];        elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
235              elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];        elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
236              elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];        elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
237              elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];        elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
238              elseif istr == 6,   Sstr = [Sstr,'(1:no),'];        elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
239              elseif istr == 7,   Sstr = [Sstr,'(1:np),'];        elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
240              else, error('Can''t handle this many dimensions!');        else, error('Can''t handle this many dimensions!');
241              end        end
242          end      end
243          eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])      eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
244          %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;      %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
245          if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end      if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
246        
247          S.attributes.(cvar)=read_att(nc{cvar});      S.attributes.(cvar)=read_att(nc{cvar});
248          % replace missing or FillValues with NaN      % replace missing or FillValues with NaN
249          attnames=fieldnames(S.attributes.(cvar));      if isfield(S.attributes.(cvar),'missing_value');
250          if ~isempty(attnames)        misval = S.attributes.(cvar).missing_value;
251            for k=1:length(attnames)        S.(cvar)(S.(cvar) == misval) = NaN;
252              if strcmp(attnames{k},'missing_value') ...      end
253                    | strcmp(attnames{k},'FillValue_')      if isfield(S.attributes.(cvar),'FillValue_');
254                S.(cvar)(S.(cvar) == S.attributes.(cvar).(attnames{k})) = NaN;        misval = S.attributes.(cvar).FillValue_;
255              end        S.(cvar)(S.(cvar) == misval) = NaN;
256            end      end
257          end      
258          end    end % for ivar
259        
260  if isempty(S)    if isempty(S)
261      error('Something didn''t work!!!');      error('Something didn''t work!!!');
262  end    end
263      
264      return

Legend:
Removed from v.1.18  
changed lines
  Added in v.1.19

  ViewVC Help
Powered by ViewVC 1.1.22