/[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.7 by enderton, Fri May 12 16:52:09 2006 UTC revision 1.19 by mlosch, Tue Nov 24 13:09:16 2009 UTC
# Line 30  function [S] = rdmnc(varargin) Line 30  function [S] = rdmnc(varargin)
30  %  Author:  Alistair Adcroft  %  Author:  Alistair Adcroft
31  %  Modifications:  Daniel Enderton  %  Modifications:  Daniel Enderton
32    
33    % $Header$
34    % $Name$
35    
36  % Initializations  % Initializations
37    dBug=0;
38  file={};  file={};
39  filepaths={};  filepaths={};
40  files={};  files={};
# Line 51  for iarg=1:nargin; Line 55  for iarg=1:nargin;
55              iters=arg;              iters=arg;
56          else          else
57              error(['The only allowed numeric argument is iterations',...              error(['The only allowed numeric argument is iterations',...
58                     ' to read in as a vector for the last arguement']);                     ' to read in as a vector for the last argument']);
59          end          end
60      end      end
61  end  end
62    if isempty(file)
63        if isempty(varlist),
64           fprintf( 'No file name in argument list\n');
65        else
66           fprintf(['No file in argument list:\n ==> ',char(varlist(1))]);
67           for i=2:size(varlist,2), fprintf([' , ',char(varlist(i))]); end
68           fprintf(' <==\n');
69        end
70        error(' check argument list !!!');
71    end
72    
73  % Create list of filenames  % Create list of filenames
74  for eachfile=file  for eachfile=file
75          filepathtemp=eachfile{:};          filepathtemp=eachfile{:};
76          indecies = find(filepathtemp=='/');          indecies = find(filepathtemp==filesep);
77          if ~isempty(indecies)          if ~isempty(indecies)
78          filepathtemp = filepathtemp(1:indecies(end));          filepathtemp = filepathtemp(1:indecies(end));
79          else          else
# Line 84  if isempty(iters) Line 98  if isempty(iters)
98          nc=netcdf(char(eachfile),'read');          nc=netcdf(char(eachfile),'read');
99          nciters = nc{'iter'}(:);          nciters = nc{'iter'}(:);
100          if isempty(nciters), nciters = nc{'T'}(:); end          if isempty(nciters), nciters = nc{'T'}(:); end
101          iters = [iters,nciters];          iters = [iters,nciters'];
102          close(nc);          close(nc);
103      end      end
104      iters = unique(iters);      iters = unique(iters');
105  end  end
106    
107  % Cycle through files for data extraction.  % Cycle through files for data extraction.
108  S.attributes=[];  S.attributes=[];
109  for ieachfile=1:length(files)  for ieachfile=1:length(files)
110      eachfile = [filepaths{ieachfile},files{ieachfile}];      eachfile = [filepaths{ieachfile},files{ieachfile}];
111        if dBug > 0, fprintf([' open: ',eachfile]); end
112      nc=netcdf(char(eachfile),'read');      nc=netcdf(char(eachfile),'read');
113      S=rdmnc_local(nc,varlist,iters,S);      S=rdmnc_local(nc,varlist,iters,S,dBug);
114      close(nc);      close(nc);
115  end  end
116    
# Line 131  function [i0,j0,fn] = findTileOffset(S); Line 146  function [i0,j0,fn] = findTileOffset(S);
146      end      end
147      %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];      %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
148    
149  function [S] = rdmnc_local(nc,varlist,iters,S)  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      [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index    if isinf(iters); iters = fiter(end); end
154      dii = dii(find(dii ~= 0));                           % Data interation index    if isnan(iters); iters = fiter; end
155          [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
156      % No variables specified? Default to all    dii = dii(find(dii ~= 0));                           % Data interation index
157      if isempty(varlist), varlist=ncnames(var(nc)); end    if dBug > 0,
158            fprintf(' ; fii='); fprintf(' %i',fii);
159      % Attributes for structure      fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
160      if iters>0; S.iters_read_from_file=iters; end    end
161      S.attributes.global=read_att(nc);      
162          % No variables specified? Default to all
163          % Read variable data    if isempty(varlist), varlist=ncnames(var(nc)); end
164          for ivar=1:size(varlist,2)    
165              % Attributes for structure
166          cvar=char(varlist{ivar});    if iters>0; S.iters_from_file=iters; end
167          if isempty(nc{cvar})    S.attributes.global=read_att(nc);
168              disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);    [pstr,netcdf_fname,ext] = fileparts(name(nc));
169              continue    if strcmp(netcdf_fname(end-3:end),'glob')
170          end      % assume it is a global file produced by gluemnc and change some
171                % attributes
172          dims = ncnames(dim(nc{cvar}));        % Dimensions      S.attributes.global.sNx = S.attributes.global.Nx;
173          adj = 0;      S.attributes.global.sNy = S.attributes.global.Ny;
174          if dims{1} == 'T'      S.attributes.global.nPx = 1;
175              if isempty(find(fii)), return, end      S.attributes.global.nSx = 1;
176                    S.attributes.global.nPy = 1;
177              tmpdata = nc{cvar}(fii,:);      S.attributes.global.nSy = 1;
178              if ismember('Zd000001' ,dims), adj = adj - 1; end      S.attributes.global.tile_number = 1;
179              if ismember('Zmd000001',dims), adj = adj - 1; end      S.attributes.global.nco_openmp_thread_number = 1;
180              it = length(dims) + adj;    end
181          else      
182              tmpdata = nc{cvar}(:);    % Read variable data
183              it = 0;    for ivar=1:size(varlist,2)
184          end      
185                cvar=char(varlist{ivar});
186          tmpdata=squeeze(permute(tmpdata,[9:-1:1]));      if isempty(nc{cvar})
187          [ni nj nk nm nn no np]=size(tmpdata);        disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
188                  continue
189          [i0,j0,fn]=findTileOffset(S);      end
190          cdim=dims{end}; if cdim(1)~='X'; i0=0; end      
191          cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end      dims = ncnames(dim(nc{cvar}));        % Dimensions
192          if length(dims)>1;      sizVar = size(nc{cvar}); nDims=length(sizVar);
193              cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end      if dims{1} == 'T'
194          else        if isempty(find(fii)), error('Iters not found'); end
195              j0=0;        it = length(dims);
196          end        tmpdata = nc{cvar}(fii,:);
197                  % leading unity dimensions get lost; add them back:
198          Sstr = '';        tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
199          for istr = 1:7      else
200              if     istr == it,  Sstr = [Sstr,'dii,'];        it = 0;
201              elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];        tmpdata = nc{cvar}(:);
202              elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];        % leading unity dimensions get lost; add them back:
203              elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];        tmpdata=reshape(tmpdata,sizVar);
204              elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];      end
205              elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];      
206              elseif istr == 6,   Sstr = [Sstr,'(1:no),'];      if dBug > 1,
207              elseif istr == 7,   Sstr = [Sstr,'(1:np),'];        fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
208              else, error('Can''t handle this many dimensions!');        fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
209              end      end
210          end      if length(dims) > 1,
211                  tmpdata=permute(tmpdata,[nDims:-1:1]);
212          eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])      end
213          %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;      if dBug > 1,
214          S.attributes.(cvar)=read_att(nc{cvar});  %         fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
215          end      end
216        [ni nj nk nm nn no np]=size(tmpdata);
217  if isempty(S)        
218        [i0,j0,fn]=findTileOffset(S);
219        cdim=dims{end}; if cdim(1)~='X'; i0=0; end
220        cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
221        if length(dims)>1;
222          cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
223        else
224          j0=0;
225        end
226        if dBug > 1,
227          fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
228        end
229        
230        Sstr = '';
231        for istr = 1:max(nDims,length(dims)),
232          if     istr == it,  Sstr = [Sstr,'dii,'];
233          elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
234          elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
235          elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
236          elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
237          elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
238          elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
239          elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
240          else, error('Can''t handle this many dimensions!');
241          end
242        end
243        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;
245        if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
246        
247        S.attributes.(cvar)=read_att(nc{cvar});
248        % replace missing or FillValues with NaN
249        if isfield(S.attributes.(cvar),'missing_value');
250          misval = S.attributes.(cvar).missing_value;
251          S.(cvar)(S.(cvar) == misval) = NaN;
252        end
253        if isfield(S.attributes.(cvar),'FillValue_');
254          misval = S.attributes.(cvar).FillValue_;
255          S.(cvar)(S.(cvar) == misval) = NaN;
256        end
257        
258      end % for ivar
259        
260      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.7  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22