utils/matlab/rdmnc.m

function [S] = rdmnc(varargin)

% Usage:
%   S=RDMNC(FILE1,FILE2,...)
%   S=RDMNC(FILE1,...,ITER)
%   S=RDMNC(FILE1,...,'VAR1','VAR2',...)
%   S=RDMNC(FILE1,...,'VAR1','VAR2',...,ITER)
%
% Input:
%   FILE1   Either a single file name (e.g. 'state.nc') or a wild-card
%           strings expanding to a group of file names (e.g. 'state.*.nc').
%           There are no assumptions about suffices (e.g. 'state.*' works).
%   VAR1    Model variable names as written in the MNC/netcdf file.
%   ITER    Vector of iterations in the MNC/netcdf files, not model time.
%
% Output:
%   S       Structure with fields corresponding to 'VAR1', 'VAR2', ...
%
% Description:
%   This function is a rudimentary wrapper for joining and reading netcdf
%   files produced by MITgcm.  It does not give the same flexibility as 
%   opening the netcdf files directly using netcdf(), but is useful for
%   quick loading of entire model fields which are distributed in multiple
%   netcdf files.
%
% Example:
%   >> S=rdmnd('state.*','XC','YC','T');
%   >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
%
%  Author:  Alistair Adcroft
%  Modifications:  Daniel Enderton

% $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.23 2011/01/27 11:42:02 mlosch Exp $
% $Name:  $

% Initializations
dBug=0;
file={};
filepaths={};
files={};
varlist={};
iters=[];

% find out if matlab's generic netcdf API is available
% if not assume that Chuck Denham's netcdf toolbox is installed
usemexcdf = isempty(which('netcdf.open'));

% Process function arguments
for iarg=1:nargin;
    arg=varargin{iarg};
    if ischar(arg)
        if isempty(dir(char(arg)))
            varlist{end+1}=arg;
        else
            file{end+1}=arg;
        end
    else
        if isempty(iters)
            iters=arg;
        else
            error(['The only allowed numeric argument is iterations',...
                   ' to read in as a vector for the last argument']);
        end
    end
end
if isempty(file)
    if isempty(varlist),
       fprintf( 'No file name in argument list\n');
    else
       fprintf(['No file in argument list:\n ==> ',char(varlist(1))]);
       for i=2:size(varlist,2), fprintf([' , ',char(varlist(i))]); end
       fprintf(' <==\n'); 
    end
    error(' check argument list !!!');
end

% Create list of filenames
for eachfile=file
        filepathtemp=eachfile{:};
        indecies = find(filepathtemp==filesep);
        if ~isempty(indecies)
        filepathtemp = filepathtemp(1:indecies(end));
        else
        filepathtemp = '';
        end
    expandedEachFile=dir(char(eachfile{1}));
    for i=1:size(expandedEachFile,1);
        if expandedEachFile(i).isdir==0
            files{end+1}=expandedEachFile(i).name;
            filepaths{end+1}=filepathtemp;
        end
    end
end


% If iterations unspecified, open all the files and make list of all the
% iterations that appear, use this iterations list for data extraction.
if isempty(iters)
    iters = [];
    for ieachfile=1:length(files)
        eachfile = [filepaths{ieachfile},files{ieachfile}];
        if usemexcdf
          nc=netcdf(char(eachfile),'read');
          nciters = nc{'iter'}(:);
          if isempty(nciters), nciters = nc{'T'}(:); end
          close(nc);
        else
          % the parser complains about "netcdf.open" when the matlab netcdf
          % API is not available, even when it is not used so we have to
          % avoid the use of "netcdf.open", etc in this function
          nciters = ncgetvar(char(eachfile),'iter');
          if isempty(nciters), nciters = ncgetvar(char(eachfile),'T'); end
        end
        iters = [iters,nciters'];
    end
    iters = unique(iters');
end

% Cycle through files for data extraction.
S.attributes=[];
for ieachfile=1:length(files)
    eachfile = [filepaths{ieachfile},files{ieachfile}];
    if dBug > 0, fprintf([' open: ',eachfile]); end
    if usemexcdf
      nc=netcdf(char(eachfile),'read');
      S=rdmnc_local(nc,varlist,iters,S,dBug);
      close(nc);
    else
      % the parser complains about "netcdf.open" when the matlab netcdf
      % API is not available, even when it is not used so we have to
      % avoid the use of "netcdf.open", etc in this function
      S=rdmnc_local_matlabAPI(char(eachfile),varlist,iters,S,dBug);
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             Local functions                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [A] = read_att(nc);
    allatt=ncnames(att(nc));
    if ~isempty(allatt)
        for attr=allatt;
            A.(char(attr))=nc.(char(attr))(:);
        end
    else
        A = 'none';
    end

function [i0,j0,fn] = findTileOffset(S);
    fn=0;
    if isfield(S,'attributes') & isfield(S.attributes,'global')
        G=S.attributes.global;
        tn=G.tile_number;
        snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
        ntx=nsx*npx;nty=nsy*npy;
        gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
        i0=snx*gbi; j0=sny*gbj;
        if isfield(S.attributes.global,'exch2_myFace')
            fn=G.exch2_myFace;
        end
    else
        i0=0;j0=0;
    end
    %[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];

function [S] = rdmnc_local(nc,varlist,iters,S,dBug)

  fiter = nc{'iter'}(:);                               % File iterations present
  if isempty(fiter), fiter = nc{'T'}(:); end
  if isinf(iters); iters = fiter(end); end
  if isnan(iters); iters = fiter; end
  [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
  dii = dii(find(dii ~= 0));                           % Data interation index
  if dBug > 0,
    fprintf(' ; fii='); fprintf(' %i',fii); 
    fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n'); 
  end
    
  % No variables specified? Default to all
  if isempty(varlist), varlist=ncnames(var(nc)); end
  
  % Attributes for structure
  if iters>0; S.iters_from_file=iters; end
  S.attributes.global=read_att(nc);
  [pstr,netcdf_fname,ext] = fileparts(name(nc));
  if strcmp(netcdf_fname(end-3:end),'glob')
    % assume it is a global file produced by gluemnc and change some
    % attributes 
    S.attributes.global.sNx = S.attributes.global.Nx;
    S.attributes.global.sNy = S.attributes.global.Ny;
    S.attributes.global.nPx = 1;
    S.attributes.global.nSx = 1;
    S.attributes.global.nPy = 1;
    S.attributes.global.nSy = 1;
    S.attributes.global.tile_number = 1;
    S.attributes.global.nco_openmp_thread_number = 1;
  end
    
  % Read variable data
  for ivar=1:size(varlist,2)
    
    cvar=char(varlist{ivar});
    if isempty(nc{cvar})
      disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
      continue
    end
    % code by Bruno Deremble: if you do not want to read all tiles these
    % modifications make the output field smaller, let us see, if it is
    % robust 
    if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
    % end code by Bruno Deremble

    dims = ncnames(dim(nc{cvar}));        % Dimensions
    sizVar = size(nc{cvar}); nDims=length(sizVar);
    if dims{1} == 'T'
      if isempty(find(fii)), error('Iters not found'); end
      it = length(dims);
      tmpdata = nc{cvar}(fii,:);
      % leading unity dimensions get lost; add them back:
      tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
    else
      it = 0;
      tmpdata = nc{cvar}(:);
      % leading unity dimensions get lost; add them back:
      tmpdata=reshape(tmpdata,sizVar);
    end
    
    if dBug > 1,
      fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
      fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it); 
    end
    if length(dims) > 1,
      tmpdata=permute(tmpdata,[nDims:-1:1]);
    end
    if dBug > 1,
%         fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
    end
    [ni nj nk nm nn no np]=size(tmpdata);
      
    [i0,j0,fn]=findTileOffset(S);
    cdim=dims{end}; if cdim(1)~='X'; i0=0; end
    cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
    if length(dims)>1;
      cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
    else
      j0=0;
    end
    if dBug > 1,
      fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
    end
    % code by Bruno Deremble: if you do not want to read all tiles these
    % modifications make the output field smaller, let us see, if it is
    % robust 
    if (firstiter)
      S.i_first.(cvar) = i0;
      S.j_first.(cvar) = j0;
    end 
    i0 = i0 - S.i_first.(cvar);
    j0 = j0 - S.j_first.(cvar);
    % end code by Bruno Deremble
    
    Sstr = '';
    for istr = 1:max(nDims,length(dims)),
      if     istr == it,  Sstr = [Sstr,'dii,'];
      elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
      elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
      elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
      elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
      elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
      elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
      elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
      else, error('Can''t handle this many dimensions!');
      end
    end
    eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
    %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
    if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
    
    S.attributes.(cvar)=read_att(nc{cvar});
    % replace missing or FillValues with NaN
    if isfield(S.attributes.(cvar),'missing_value');
      misval = S.attributes.(cvar).missing_value;
      S.(cvar)(S.(cvar) == misval) = NaN;
    end
    if isfield(S.attributes.(cvar),'FillValue_');
      misval = S.attributes.(cvar).FillValue_;
      S.(cvar)(S.(cvar) == misval) = NaN;
    end
    
  end % for ivar
    
  if isempty(S)
    error('Something didn''t work!!!');
  end
  
  return

function [S] = rdmnc_local_matlabAPI(fname,varlist,iters,S,dBug)

  fiter = ncgetvar(fname,'iter');                     % File iterations present
  if isempty(fiter), fiter = ncgetvar(fname,'T'); end
  if isinf(iters); iters = fiter(end); end
  if isnan(iters); iters = fiter; end
  [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
  dii = dii(find(dii ~= 0));                           % Data interation index
  if dBug > 0,
    fprintf(' ; fii='); fprintf(' %i',fii); 
    fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n'); 
  end

  % now open the file for reading
  nc = netcdf.open(fname,'NC_NOWRITE');
  % get basic information about netcdf file
  [ndims nvars natts recdim] = netcdf.inq(nc);

  % No variables specified? Default to all
  if isempty(varlist), 
    for k=0:nvars-1
      varlist{k+1} = netcdf.inqVar(nc,k);
    end
  end
  
  % Attributes for structure
  if iters>0; S.iters_from_file=iters; end
  S.attributes.global=ncgetatt(nc,'global');
  [pstr,netcdf_fname,ext] = fileparts(fname);
  if strcmp(netcdf_fname(end-3:end),'glob')
    % assume it is a global file produced by gluemnc and change some
    % attributes 
    S.attributes.global.sNx = S.attributes.global.Nx;
    S.attributes.global.sNy = S.attributes.global.Ny;
    S.attributes.global.nPx = 1;
    S.attributes.global.nSx = 1;
    S.attributes.global.nPy = 1;
    S.attributes.global.nSy = 1;
    S.attributes.global.tile_number = 1;
    S.attributes.global.nco_openmp_thread_number = 1;
  end
    
  % Read variable data
  for ivar=1:size(varlist,2)
    
    cvar=char(varlist{ivar});
    varid=ncfindvarid(nc,cvar);
    if isempty(varid)
      disp(['No such variable ''',cvar,''' in MNC file ',fname]);
      continue
    end
    % code by Bruno Deremble: if you do not want to read all tiles these
    % modifications make the output field smaller, let us see, if it is
    % robust 
    if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
    % end code by Bruno Deremble
    
    [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
    % does this variable contain a record (unlimited) dimension?
    [isrecvar,recpos] = ismember(recdim,dimids);
    
    % Dimensions
    clear sizVar dims
    for k=1:length(dimids)
      [dims{k}, sizVar(k)] = netcdf.inqDim(nc,dimids(k));
    end
    nDims=length(sizVar);
    if isrecvar
      if isempty(find(fii)), error('Iters not found'); end
      it = length(dims);
      if length(dimids) == 1
        % just a time or iteration variable, this will result in a vector
        % and requires special treatment
        icount=1;
        tmpdata = zeros(length(fii),1);
        for k=1:length(fii)
          istart = fii(k)-1;
          tmpdata(k) = netcdf.getVar(nc,varid,istart,icount,'double');
        end
      else
        % from now on we assume that the record dimension is always the
        % last dimension. This may not always be the case
        if recpos ~= nDims
          error(sprintf('%s\n%s\n%s%s%i%s%i', ...
                        ['The current code assumes that the record ' ...
                         'dimension is the last dimension,'], ...
                        'this is not the case for variable', cvar, ...
                        ': nDims = ', nDims,  ...
                        ', position of recDim = ', recpos))
        end
        istart = zeros(1,it); % indexing starts a 0
        icount = sizVar;
        % we always want to get only on time slice at a time 
        icount(recpos) = 1; 
        % make your life simpler by putting the time dimension first
        tmpdata = zeros([length(fii) sizVar(1:end-1)]);
        for k=1:length(fii)
          istart(recpos) = fii(k)-1; % indexing starts at 0
          tmp = netcdf.getVar(nc,varid,istart,icount,'double');
          tmpdata(k,:) = tmp(:);
        end
        % move time dimension to the end ...
        tmpdata = shiftdim(tmpdata,1);
        % ... and restore original shape
        tmpdata = reshape(tmpdata,[sizVar(1:end-1) length(fii)]);
      end
    else
      it = 0;
      tmpdata = netcdf.getVar(nc,varid,'double');
    end
    %
    if dBug > 1,
      fprintf(['  var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',sizVar);
      fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it); 
    end
    [ni nj nk nm nn no np]=size(tmpdata);
    %
    [i0,j0,fn]=findTileOffset(S);
    cdim=dims{1}; if cdim(1)~='X'; i0=0; end
    cdim=dims{1}; if cdim(1)=='Y'; i0=j0; j0=0; end
    if length(dims)>1;
      cdim=dims{2}; if cdim(1)~='Y'; j0=0; end
    else
      j0=0;
    end
    if dBug > 1,
      fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
    end
    % code by Bruno Deremble: if you do not want to read all tiles these
    % modifications make the output field smaller, let us see, if it is
    % robust 
    if (firstiter)
      S.i_first.(cvar) = i0;
      S.j_first.(cvar) = j0;
    end 
    i0 = i0 - S.i_first.(cvar);
    j0 = j0 - S.j_first.(cvar);
    % end code by Bruno Deremble

    Sstr = '';
    for istr = 1:max(nDims,length(dims)),
      if     istr == it,  Sstr = [Sstr,'dii,'];
      elseif istr == 1,   Sstr = [Sstr,'i0+(1:ni),'];
      elseif istr == 2,   Sstr = [Sstr,'j0+(1:nj),'];
      elseif istr == 3,   Sstr = [Sstr,'(1:nk),'];
      elseif istr == 4,   Sstr = [Sstr,'(1:nm),'];
      elseif istr == 5,   Sstr = [Sstr,'(1:nn),'];
      elseif istr == 6,   Sstr = [Sstr,'(1:no),'];
      elseif istr == 7,   Sstr = [Sstr,'(1:np),'];
      else, error('Can''t handle this many dimensions!');
      end
    end
    eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
    %S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
    if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
    %
    S.attributes.(cvar)=ncgetatt(nc,cvar);
    % replace missing or FillValues with NaN
    if isfield(S.attributes.(cvar),'missing_value');
      misval = S.attributes.(cvar).missing_value;
      S.(cvar)(S.(cvar) == misval) = NaN;
    end
    if isfield(S.attributes.(cvar),'FillValue_');
      misval = S.attributes.(cvar).FillValue_;
      S.(cvar)(S.(cvar) == misval) = NaN;
    end
    
  end % for ivar

  % close the file
  netcdf.close(nc);
  
  if isempty(S)
    error('Something didn''t work!!!');
  end
  
  return

function vf = ncgetvar(fname,varname)
% read a netcdf variable
  
  nc=netcdf.open(fname,'NC_NOWRITE');
  % find out basics about the files
  [ndims nvars natts dimm] = netcdf.inq(nc);
  vf = [];
  varid = [];
  for k=0:nvars-1
    if strcmp(netcdf.inqVar(nc,k),varname)
      varid = netcdf.inqVarID(nc,varname);
    end
  end
  if ~isempty(varid); 
    [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
    % get data
    vf = double(netcdf.getVar(nc,varid));
  else
    % do nothing
  end
  netcdf.close(nc);
  
  return

function misval = ncgetmisval(nc,varid)

  [varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
  misval = [];
  for k=0:natts-1
    attn = netcdf.inqAttName(nc,varid,k);
    if strcmp(attn,'missing_value') | strcmp(attn,'_FillValue')
      misval = double(netcdf.getAtt(nc,varid,attname));
    end
  end
  
function A = ncgetatt(nc,varname)
% get all attributes and put them into a struct
  
% 1. get global properties of file
  [ndims nvars natts dimm] = netcdf.inq(nc);

  % get variable ID and properties
  if strcmp('global',varname)
    % "variable" is global
    varid = netcdf.getConstant('NC_GLOBAL');
  else
    % find variable ID and properties
    varid = ncfindvarid(nc,varname);
    if ~isempty(varid)
      [varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
    else
      warning(sprintf('variable %s not found',varname))
    end
  end

  if natts > 1
    for k=0:natts-1
      attn = netcdf.inqAttName(nc,varid,k);
      [xtype attlen]=netcdf.inqAtt(nc,varid,attn);
      attval = netcdf.getAtt(nc,varid,attn);
      if ~ischar(attval)
        attval = double(attval);
      end
      if strcmp(attn,'_FillValue')
        % matlab does not allow variable names to begin with an
        % underscore ("_"), so we have to do change the name of this
        % obsolete attribute.
        A.FillValue_=attval;
      else
        A.(char(attn))=attval;
      end
    end
  else
      A = 'none';
  end
  
  return

  
function varid = ncfindvarid(nc,varname)

  [ndims nvars natts dimm] = netcdf.inq(nc);
  varid=[];
  for k=0:nvars-1
    if strcmp(netcdf.inqVar(nc,k),varname);
      varid = netcdf.inqVarID(nc,varname);
      break
    end
  end

  return
1	enderton	1.6	function [S] = rdmnc(varargin)
2	enderton	1.3
3			% Usage:
4			% S=RDMNC(FILE1,FILE2,...)
5			% S=RDMNC(FILE1,...,ITER)
6			% S=RDMNC(FILE1,...,'VAR1','VAR2',...)
7			% S=RDMNC(FILE1,...,'VAR1','VAR2',...,ITER)
8			%
9			% Input:
10			% FILE1 Either a single file name (e.g. 'state.nc') or a wild-card
11			% strings expanding to a group of file names (e.g. 'state.*.nc').
12			% There are no assumptions about suffices (e.g. 'state.*' works).
13			% VAR1 Model variable names as written in the MNC/netcdf file.
14			% ITER Vector of iterations in the MNC/netcdf files, not model time.
15	adcroft	1.1	%
16	enderton	1.3	% Output:
17			% S Structure with fields corresponding to 'VAR1', 'VAR2', ...
18	adcroft	1.1	%
19	enderton	1.3	% Description:
20			% This function is a rudimentary wrapper for joining and reading netcdf
21			% files produced by MITgcm. It does not give the same flexibility as
22			% opening the netcdf files directly using netcdf(), but is useful for
23			% quick loading of entire model fields which are distributed in multiple
24			% netcdf files.
25	adcroft	1.1	%
26	enderton	1.3	% Example:
27			% >> S=rdmnd('state.*','XC','YC','T');
28			% >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
29	adcroft	1.1	%
30	enderton	1.4	% Author: Alistair Adcroft
31	enderton	1.3	% Modifications: Daniel Enderton
32	adcroft	1.1
33	mlosch	1.24	% $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.23 2011/01/27 11:42:02 mlosch Exp $
34	jmc	1.12	% $Name: $
35
36	enderton	1.3	% Initializations
37	jmc	1.13	dBug=0;
38	adcroft	1.1	file={};
39	enderton	1.3	filepaths={};
40	adcroft	1.1	files={};
41			varlist={};
42	enderton	1.3	iters=[];
43	adcroft	1.1
44	mlosch	1.20	% find out if matlab's generic netcdf API is available
45			% if not assume that Chuck Denham's netcdf toolbox is installed
46			usemexcdf = isempty(which('netcdf.open'));
47
48	adcroft	1.1	% Process function arguments
49			for iarg=1:nargin;
50	enderton	1.3	arg=varargin{iarg};
51			if ischar(arg)
52			if isempty(dir(char(arg)))
53			varlist{end+1}=arg;
54			else
55			file{end+1}=arg;
56			end
57			else
58			if isempty(iters)
59			iters=arg;
60			else
61			error(['The only allowed numeric argument is iterations',...
62	jmc	1.14	' to read in as a vector for the last argument']);
63	enderton	1.3	end
64			end
65	adcroft	1.1	end
66	jmc	1.14	if isempty(file)
67			if isempty(varlist),
68			fprintf( 'No file name in argument list\n');
69			else
70			fprintf(['No file in argument list:\n ==> ',char(varlist(1))]);
71			for i=2:size(varlist,2), fprintf([' , ',char(varlist(i))]); end
72			fprintf(' <==\n');
73			end
74			error(' check argument list !!!');
75			end
76	adcroft	1.1
77			% Create list of filenames
78			for eachfile=file
79	enderton	1.3	filepathtemp=eachfile{:};
80	mlosch	1.17	indecies = find(filepathtemp==filesep);
81	enderton	1.3	if ~isempty(indecies)
82			filepathtemp = filepathtemp(1:indecies(end));
83			else
84			filepathtemp = '';
85			end
86			expandedEachFile=dir(char(eachfile{1}));
87			for i=1:size(expandedEachFile,1);
88			if expandedEachFile(i).isdir==0
89			files{end+1}=expandedEachFile(i).name;
90			filepaths{end+1}=filepathtemp;
91			end
92			end
93	adcroft	1.1	end
94
95
96	enderton	1.3	% If iterations unspecified, open all the files and make list of all the
97			% iterations that appear, use this iterations list for data extraction.
98			if isempty(iters)
99			iters = [];
100			for ieachfile=1:length(files)
101			eachfile = [filepaths{ieachfile},files{ieachfile}];
102	mlosch	1.20	if usemexcdf
103			nc=netcdf(char(eachfile),'read');
104			nciters = nc{'iter'}(:);
105			if isempty(nciters), nciters = nc{'T'}(:); end
106	mlosch	1.21	close(nc);
107	mlosch	1.20	else
108	mlosch	1.21	% the parser complains about "netcdf.open" when the matlab netcdf
109			% API is not available, even when it is not used so we have to
110			% avoid the use of "netcdf.open", etc in this function
111			nciters = ncgetvar(char(eachfile),'iter');
112			if isempty(nciters), nciters = ncgetvar(char(eachfile),'T'); end
113	mlosch	1.20	end
114	jmc	1.15	iters = [iters,nciters'];
115	enderton	1.3	end
116	jmc	1.15	iters = unique(iters');
117	adcroft	1.1	end
118
119	enderton	1.3	% Cycle through files for data extraction.
120			S.attributes=[];
121			for ieachfile=1:length(files)
122			eachfile = [filepaths{ieachfile},files{ieachfile}];
123	jmc	1.13	if dBug > 0, fprintf([' open: ',eachfile]); end
124	mlosch	1.20	if usemexcdf
125			nc=netcdf(char(eachfile),'read');
126			S=rdmnc_local(nc,varlist,iters,S,dBug);
127			close(nc);
128			else
129	mlosch	1.21	% the parser complains about "netcdf.open" when the matlab netcdf
130			% API is not available, even when it is not used so we have to
131			% avoid the use of "netcdf.open", etc in this function
132			S=rdmnc_local_matlabAPI(char(eachfile),varlist,iters,S,dBug);
133	mlosch	1.20	end
134	adcroft	1.1	end
135
136	mlosch	1.20	return
137	enderton	1.5
138	enderton	1.3	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139			% Local functions %
140			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141	adcroft	1.1
142	enderton	1.3	function [A] = read_att(nc);
143			allatt=ncnames(att(nc));
144			if ~isempty(allatt)
145			for attr=allatt;
146			A.(char(attr))=nc.(char(attr))(:);
147			end
148			else
149			A = 'none';
150			end
151	adcroft	1.1
152	enderton	1.3	function [i0,j0,fn] = findTileOffset(S);
153			fn=0;
154			if isfield(S,'attributes') & isfield(S.attributes,'global')
155			G=S.attributes.global;
156			tn=G.tile_number;
157			snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
158			ntx=nsxnpx;nty=nsynpy;
159			gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
160			i0=snxgbi; j0=snygbj;
161			if isfield(S.attributes.global,'exch2_myFace')
162			fn=G.exch2_myFace;
163			end
164			else
165			i0=0;j0=0;
166			end
167			%[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
168
169	jmc	1.13	function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
170	enderton	1.3
171	mlosch	1.19	fiter = nc{'iter'}(:); % File iterations present
172			if isempty(fiter), fiter = nc{'T'}(:); end
173			if isinf(iters); iters = fiter(end); end
174			if isnan(iters); iters = fiter; end
175			[fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
176			dii = dii(find(dii ~= 0)); % Data interation index
177			if dBug > 0,
178			fprintf(' ; fii='); fprintf(' %i',fii);
179			fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
180			end
181
182			% No variables specified? Default to all
183			if isempty(varlist), varlist=ncnames(var(nc)); end
184
185			% Attributes for structure
186			if iters>0; S.iters_from_file=iters; end
187			S.attributes.global=read_att(nc);
188			[pstr,netcdf_fname,ext] = fileparts(name(nc));
189			if strcmp(netcdf_fname(end-3:end),'glob')
190			% assume it is a global file produced by gluemnc and change some
191			% attributes
192			S.attributes.global.sNx = S.attributes.global.Nx;
193			S.attributes.global.sNy = S.attributes.global.Ny;
194			S.attributes.global.nPx = 1;
195			S.attributes.global.nSx = 1;
196			S.attributes.global.nPy = 1;
197			S.attributes.global.nSy = 1;
198			S.attributes.global.tile_number = 1;
199			S.attributes.global.nco_openmp_thread_number = 1;
200			end
201
202			% Read variable data
203			for ivar=1:size(varlist,2)
204
205			cvar=char(varlist{ivar});
206			if isempty(nc{cvar})
207			disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
208			continue
209			end
210	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
211			% modifications make the output field smaller, let us see, if it is
212			% robust
213			if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
214			% end code by Bruno Deremble
215
216	mlosch	1.19	dims = ncnames(dim(nc{cvar})); % Dimensions
217			sizVar = size(nc{cvar}); nDims=length(sizVar);
218			if dims{1} == 'T'
219			if isempty(find(fii)), error('Iters not found'); end
220			it = length(dims);
221			tmpdata = nc{cvar}(fii,:);
222			% leading unity dimensions get lost; add them back:
223			tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
224			else
225			it = 0;
226			tmpdata = nc{cvar}(:);
227			% leading unity dimensions get lost; add them back:
228			tmpdata=reshape(tmpdata,sizVar);
229			end
230
231			if dBug > 1,
232			fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
233			fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
234			end
235			if length(dims) > 1,
236			tmpdata=permute(tmpdata,[nDims:-1:1]);
237			end
238			if dBug > 1,
239	jmc	1.13	% fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
240	mlosch	1.19	end
241			[ni nj nk nm nn no np]=size(tmpdata);
242
243			[i0,j0,fn]=findTileOffset(S);
244			cdim=dims{end}; if cdim(1)~='X'; i0=0; end
245			cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
246			if length(dims)>1;
247			cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
248			else
249			j0=0;
250			end
251			if dBug > 1,
252			fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
253			end
254	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
255			% modifications make the output field smaller, let us see, if it is
256			% robust
257			if (firstiter)
258			S.i_first.(cvar) = i0;
259			S.j_first.(cvar) = j0;
260			end
261			i0 = i0 - S.i_first.(cvar);
262			j0 = j0 - S.j_first.(cvar);
263			% end code by Bruno Deremble
264	mlosch	1.19
265			Sstr = '';
266			for istr = 1:max(nDims,length(dims)),
267			if istr == it, Sstr = [Sstr,'dii,'];
268			elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
269			elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
270			elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
271			elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
272			elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
273			elseif istr == 6, Sstr = [Sstr,'(1:no),'];
274			elseif istr == 7, Sstr = [Sstr,'(1:np),'];
275			else, error('Can''t handle this many dimensions!');
276			end
277			end
278			eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
279			%S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
280			if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
281
282			S.attributes.(cvar)=read_att(nc{cvar});
283			% replace missing or FillValues with NaN
284			if isfield(S.attributes.(cvar),'missing_value');
285			misval = S.attributes.(cvar).missing_value;
286			S.(cvar)(S.(cvar) == misval) = NaN;
287			end
288			if isfield(S.attributes.(cvar),'FillValue_');
289			misval = S.attributes.(cvar).FillValue_;
290			S.(cvar)(S.(cvar) == misval) = NaN;
291			end
292
293			end % for ivar
294
295			if isempty(S)
296	enderton	1.3	error('Something didn''t work!!!');
297	mlosch	1.19	end
298
299			return
300	mlosch	1.20
301	mlosch	1.21	function [S] = rdmnc_local_matlabAPI(fname,varlist,iters,S,dBug)
302	mlosch	1.20
303	mlosch	1.21	fiter = ncgetvar(fname,'iter'); % File iterations present
304			if isempty(fiter), fiter = ncgetvar(fname,'T'); end
305	mlosch	1.20	if isinf(iters); iters = fiter(end); end
306			if isnan(iters); iters = fiter; end
307			[fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
308			dii = dii(find(dii ~= 0)); % Data interation index
309			if dBug > 0,
310			fprintf(' ; fii='); fprintf(' %i',fii);
311			fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
312			end
313
314	mlosch	1.21	% now open the file for reading
315			nc = netcdf.open(fname,'NC_NOWRITE');
316	mlosch	1.20	% get basic information about netcdf file
317	mlosch	1.22	[ndims nvars natts recdim] = netcdf.inq(nc);
318	mlosch	1.20
319			% No variables specified? Default to all
320			if isempty(varlist),
321			for k=0:nvars-1
322			varlist{k+1} = netcdf.inqVar(nc,k);
323			end
324			end
325
326			% Attributes for structure
327			if iters>0; S.iters_from_file=iters; end
328			S.attributes.global=ncgetatt(nc,'global');
329			[pstr,netcdf_fname,ext] = fileparts(fname);
330			if strcmp(netcdf_fname(end-3:end),'glob')
331			% assume it is a global file produced by gluemnc and change some
332			% attributes
333			S.attributes.global.sNx = S.attributes.global.Nx;
334			S.attributes.global.sNy = S.attributes.global.Ny;
335			S.attributes.global.nPx = 1;
336			S.attributes.global.nSx = 1;
337			S.attributes.global.nPy = 1;
338			S.attributes.global.nSy = 1;
339			S.attributes.global.tile_number = 1;
340			S.attributes.global.nco_openmp_thread_number = 1;
341			end
342
343			% Read variable data
344			for ivar=1:size(varlist,2)
345
346			cvar=char(varlist{ivar});
347			varid=ncfindvarid(nc,cvar);
348			if isempty(varid)
349			disp(['No such variable ''',cvar,''' in MNC file ',fname]);
350			continue
351			end
352	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
353			% modifications make the output field smaller, let us see, if it is
354			% robust
355			if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
356			% end code by Bruno Deremble
357	mlosch	1.20
358			[varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
359	mlosch	1.22	% does this variable contain a record (unlimited) dimension?
360			[isrecvar,recpos] = ismember(recdim,dimids);
361	mlosch	1.20
362			% Dimensions
363			clear sizVar dims
364			for k=1:length(dimids)
365			[dims{k}, sizVar(k)] = netcdf.inqDim(nc,dimids(k));
366			end
367			nDims=length(sizVar);
368	mlosch	1.22	if isrecvar
369	mlosch	1.20	if isempty(find(fii)), error('Iters not found'); end
370			it = length(dims);
371	mlosch	1.22	if length(dimids) == 1
372			% just a time or iteration variable, this will result in a vector
373			% and requires special treatment
374			icount=1;
375			tmpdata = zeros(length(fii),1);
376			for k=1:length(fii)
377			istart = fii(k)-1;
378			tmpdata(k) = netcdf.getVar(nc,varid,istart,icount,'double');
379			end
380			else
381			% from now on we assume that the record dimension is always the
382			% last dimension. This may not always be the case
383			if recpos ~= nDims
384			error(sprintf('%s\n%s\n%s%s%i%s%i', ...
385			['The current code assumes that the record ' ...
386			'dimension is the last dimension,'], ...
387			'this is not the case for variable', cvar, ...
388			': nDims = ', nDims, ...
389			', position of recDim = ', recpos))
390			end
391			istart = zeros(1,it); % indexing starts a 0
392			icount = sizVar;
393			% we always want to get only on time slice at a time
394			icount(recpos) = 1;
395			% make your life simpler by putting the time dimension first
396			tmpdata = zeros([length(fii) sizVar(1:end-1)]);
397			for k=1:length(fii)
398			istart(recpos) = fii(k)-1; % indexing starts at 0
399			tmp = netcdf.getVar(nc,varid,istart,icount,'double');
400			tmpdata(k,:) = tmp(:);
401			end
402			% move time dimension to the end ...
403			tmpdata = shiftdim(tmpdata,1);
404			% ... and restore original shape
405			tmpdata = reshape(tmpdata,[sizVar(1:end-1) length(fii)]);
406			end
407	mlosch	1.20	else
408			it = 0;
409			tmpdata = netcdf.getVar(nc,varid,'double');
410			end
411	mlosch	1.22	%
412	mlosch	1.20	if dBug > 1,
413			fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',sizVar);
414			fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
415			end
416			[ni nj nk nm nn no np]=size(tmpdata);
417	mlosch	1.22	%
418	mlosch	1.20	[i0,j0,fn]=findTileOffset(S);
419			cdim=dims{1}; if cdim(1)~='X'; i0=0; end
420			cdim=dims{1}; if cdim(1)=='Y'; i0=j0; j0=0; end
421			if length(dims)>1;
422			cdim=dims{2}; if cdim(1)~='Y'; j0=0; end
423			else
424			j0=0;
425			end
426			if dBug > 1,
427			fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
428			end
429	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
430			% modifications make the output field smaller, let us see, if it is
431			% robust
432			if (firstiter)
433			S.i_first.(cvar) = i0;
434			S.j_first.(cvar) = j0;
435			end
436			i0 = i0 - S.i_first.(cvar);
437			j0 = j0 - S.j_first.(cvar);
438			% end code by Bruno Deremble
439
440	mlosch	1.20	Sstr = '';
441			for istr = 1:max(nDims,length(dims)),
442			if istr == it, Sstr = [Sstr,'dii,'];
443			elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
444			elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
445			elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
446			elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
447			elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
448			elseif istr == 6, Sstr = [Sstr,'(1:no),'];
449			elseif istr == 7, Sstr = [Sstr,'(1:np),'];
450			else, error('Can''t handle this many dimensions!');
451			end
452			end
453			eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
454			%S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
455			if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
456	mlosch	1.22	%
457	mlosch	1.20	S.attributes.(cvar)=ncgetatt(nc,cvar);
458			% replace missing or FillValues with NaN
459			if isfield(S.attributes.(cvar),'missing_value');
460			misval = S.attributes.(cvar).missing_value;
461			S.(cvar)(S.(cvar) == misval) = NaN;
462			end
463			if isfield(S.attributes.(cvar),'FillValue_');
464			misval = S.attributes.(cvar).FillValue_;
465			S.(cvar)(S.(cvar) == misval) = NaN;
466			end
467
468			end % for ivar
469	mlosch	1.21
470			% close the file
471			netcdf.close(nc);
472
473	mlosch	1.20	if isempty(S)
474			error('Something didn''t work!!!');
475			end
476
477			return
478
479	mlosch	1.21	function vf = ncgetvar(fname,varname)
480	mlosch	1.20	% read a netcdf variable
481
482	mlosch	1.21	nc=netcdf.open(fname,'NC_NOWRITE');
483			% find out basics about the files
484	mlosch	1.20	[ndims nvars natts dimm] = netcdf.inq(nc);
485			vf = [];
486			varid = [];
487			for k=0:nvars-1
488			if strcmp(netcdf.inqVar(nc,k),varname)
489	mlosch	1.24	varid = netcdf.inqVarID(nc,varname);
490	mlosch	1.20	end
491			end
492			if ~isempty(varid);
493			[varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
494			% get data
495			vf = double(netcdf.getVar(nc,varid));
496			else
497			% do nothing
498			end
499	mlosch	1.21	netcdf.close(nc);
500
501	mlosch	1.20	return
502
503			function misval = ncgetmisval(nc,varid)
504
505			[varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
506			misval = [];
507			for k=0:natts-1
508			attn = netcdf.inqAttName(nc,varid,k);
509			if strcmp(attn,'missing_value') \| strcmp(attn,'_FillValue')
510			misval = double(netcdf.getAtt(nc,varid,attname));
511			end
512			end
513
514			function A = ncgetatt(nc,varname)
515			% get all attributes and put them into a struct
516
517			% 1. get global properties of file
518			[ndims nvars natts dimm] = netcdf.inq(nc);
519
520			% get variable ID and properties
521			if strcmp('global',varname)
522			% "variable" is global
523			varid = netcdf.getConstant('NC_GLOBAL');
524			else
525			% find variable ID and properties
526			varid = ncfindvarid(nc,varname);
527			if ~isempty(varid)
528			[varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
529			else
530			warning(sprintf('variable %s not found',varname))
531			end
532			end
533
534			if natts > 1
535			for k=0:natts-1
536			attn = netcdf.inqAttName(nc,varid,k);
537			[xtype attlen]=netcdf.inqAtt(nc,varid,attn);
538			attval = netcdf.getAtt(nc,varid,attn);
539			if ~ischar(attval)
540			attval = double(attval);
541			end
542	mlosch	1.23	if strcmp(attn,'_FillValue')
543			% matlab does not allow variable names to begin with an
544			% underscore ("_"), so we have to do change the name of this
545			% obsolete attribute.
546			A.FillValue_=attval;
547			else
548			A.(char(attn))=attval;
549			end
550	mlosch	1.20	end
551			else
552			A = 'none';
553			end
554
555			return
556
557
558			function varid = ncfindvarid(nc,varname)
559
560			[ndims nvars natts dimm] = netcdf.inq(nc);
561			varid=[];
562			for k=0:nvars-1
563			if strcmp(netcdf.inqVar(nc,k),varname);
564	mlosch	1.24	varid = netcdf.inqVarID(nc,varname);
565	mlosch	1.20	break
566			end
567			end
568
569			return