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.26 2011/06/28 09:36:42 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;
            i0=G.exch2_txGlobalo -1; j0=G.exch2_tyGlobalo -1;
        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.attributes.i_first.(cvar) = i0;
      S.attributes.j_first.(cvar) = j0;
    end
    i0 = i0 - S.attributes.i_first.(cvar);
    j0 = j0 - S.attributes.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.attributes.i_first.(cvar) = i0;
      S.attributes.j_first.(cvar) = j0;
    end
    i0 = i0 - S.attributes.i_first.(cvar);
    j0 = j0 - S.attributes.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	jmc	1.27	% files produced by MITgcm. It does not give the same flexibility as
22	enderton	1.3	% 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	jmc	1.27	% $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.26 2011/06/28 09:36:42 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	jmc	1.27	fprintf(' <==\n');
73	jmc	1.14	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	jmc	1.27	i0=G.exch2_txGlobalo -1; j0=G.exch2_tyGlobalo -1;
164	enderton	1.3	end
165			else
166			i0=0;j0=0;
167			end
168			%[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
169
170	jmc	1.13	function [S] = rdmnc_local(nc,varlist,iters,S,dBug)
171	enderton	1.3
172	mlosch	1.19	fiter = nc{'iter'}(:); % File iterations present
173			if isempty(fiter), fiter = nc{'T'}(:); end
174			if isinf(iters); iters = fiter(end); end
175			if isnan(iters); iters = fiter; end
176			[fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
177			dii = dii(find(dii ~= 0)); % Data interation index
178			if dBug > 0,
179	jmc	1.27	fprintf(' ; fii='); fprintf(' %i',fii);
180			fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
181	mlosch	1.19	end
182	jmc	1.27
183	mlosch	1.19	% No variables specified? Default to all
184			if isempty(varlist), varlist=ncnames(var(nc)); end
185	jmc	1.27
186	mlosch	1.19	% Attributes for structure
187			if iters>0; S.iters_from_file=iters; end
188			S.attributes.global=read_att(nc);
189			[pstr,netcdf_fname,ext] = fileparts(name(nc));
190			if strcmp(netcdf_fname(end-3:end),'glob')
191			% assume it is a global file produced by gluemnc and change some
192	jmc	1.27	% attributes
193	mlosch	1.19	S.attributes.global.sNx = S.attributes.global.Nx;
194			S.attributes.global.sNy = S.attributes.global.Ny;
195			S.attributes.global.nPx = 1;
196			S.attributes.global.nSx = 1;
197			S.attributes.global.nPy = 1;
198			S.attributes.global.nSy = 1;
199			S.attributes.global.tile_number = 1;
200			S.attributes.global.nco_openmp_thread_number = 1;
201			end
202	jmc	1.27
203	mlosch	1.19	% Read variable data
204			for ivar=1:size(varlist,2)
205	jmc	1.27
206	mlosch	1.19	cvar=char(varlist{ivar});
207			if isempty(nc{cvar})
208			disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
209			continue
210			end
211	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
212			% modifications make the output field smaller, let us see, if it is
213	jmc	1.27	% robust
214	mlosch	1.25	if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
215			% end code by Bruno Deremble
216
217	mlosch	1.19	dims = ncnames(dim(nc{cvar})); % Dimensions
218			sizVar = size(nc{cvar}); nDims=length(sizVar);
219			if dims{1} == 'T'
220			if isempty(find(fii)), error('Iters not found'); end
221			it = length(dims);
222			tmpdata = nc{cvar}(fii,:);
223			% leading unity dimensions get lost; add them back:
224			tmpdata=reshape(tmpdata,[length(fii) sizVar(2:end)]);
225			else
226			it = 0;
227			tmpdata = nc{cvar}(:);
228			% leading unity dimensions get lost; add them back:
229			tmpdata=reshape(tmpdata,sizVar);
230			end
231	jmc	1.27
232	mlosch	1.19	if dBug > 1,
233			fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',size(nc{cvar}));
234	jmc	1.27	fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
235	mlosch	1.19	end
236			if length(dims) > 1,
237			tmpdata=permute(tmpdata,[nDims:-1:1]);
238			end
239			if dBug > 1,
240	jmc	1.13	% fprintf('(tmpdata:');fprintf(' %i',size(tmpdata)); fprintf(')');
241	mlosch	1.19	end
242			[ni nj nk nm nn no np]=size(tmpdata);
243	jmc	1.27
244	mlosch	1.19	[i0,j0,fn]=findTileOffset(S);
245			cdim=dims{end}; if cdim(1)~='X'; i0=0; end
246			cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
247			if length(dims)>1;
248			cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
249			else
250			j0=0;
251			end
252			if dBug > 1,
253			fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
254			end
255	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
256			% modifications make the output field smaller, let us see, if it is
257	jmc	1.27	% robust
258	mlosch	1.25	if (firstiter)
259	mlosch	1.26	S.attributes.i_first.(cvar) = i0;
260			S.attributes.j_first.(cvar) = j0;
261	jmc	1.27	end
262	mlosch	1.26	i0 = i0 - S.attributes.i_first.(cvar);
263			j0 = j0 - S.attributes.j_first.(cvar);
264	mlosch	1.25	% end code by Bruno Deremble
265	jmc	1.27
266	mlosch	1.19	Sstr = '';
267			for istr = 1:max(nDims,length(dims)),
268			if istr == it, Sstr = [Sstr,'dii,'];
269			elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
270			elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
271			elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
272			elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
273			elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
274			elseif istr == 6, Sstr = [Sstr,'(1:no),'];
275			elseif istr == 7, Sstr = [Sstr,'(1:np),'];
276			else, error('Can''t handle this many dimensions!');
277			end
278			end
279			eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
280			%S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
281			if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
282	jmc	1.27
283	mlosch	1.19	S.attributes.(cvar)=read_att(nc{cvar});
284			% replace missing or FillValues with NaN
285			if isfield(S.attributes.(cvar),'missing_value');
286			misval = S.attributes.(cvar).missing_value;
287			S.(cvar)(S.(cvar) == misval) = NaN;
288			end
289			if isfield(S.attributes.(cvar),'FillValue_');
290			misval = S.attributes.(cvar).FillValue_;
291			S.(cvar)(S.(cvar) == misval) = NaN;
292			end
293	jmc	1.27
294	mlosch	1.19	end % for ivar
295	jmc	1.27
296	mlosch	1.19	if isempty(S)
297	enderton	1.3	error('Something didn''t work!!!');
298	mlosch	1.19	end
299	jmc	1.27
300	mlosch	1.19	return
301	mlosch	1.20
302	mlosch	1.21	function [S] = rdmnc_local_matlabAPI(fname,varlist,iters,S,dBug)
303	mlosch	1.20
304	mlosch	1.21	fiter = ncgetvar(fname,'iter'); % File iterations present
305			if isempty(fiter), fiter = ncgetvar(fname,'T'); end
306	mlosch	1.20	if isinf(iters); iters = fiter(end); end
307			if isnan(iters); iters = fiter; end
308			[fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
309			dii = dii(find(dii ~= 0)); % Data interation index
310			if dBug > 0,
311	jmc	1.27	fprintf(' ; fii='); fprintf(' %i',fii);
312			fprintf(' ; dii='); fprintf(' %i',dii); fprintf(' \n');
313	mlosch	1.20	end
314
315	mlosch	1.21	% now open the file for reading
316			nc = netcdf.open(fname,'NC_NOWRITE');
317	mlosch	1.20	% get basic information about netcdf file
318	mlosch	1.22	[ndims nvars natts recdim] = netcdf.inq(nc);
319	mlosch	1.20
320			% No variables specified? Default to all
321	jmc	1.27	if isempty(varlist),
322	mlosch	1.20	for k=0:nvars-1
323			varlist{k+1} = netcdf.inqVar(nc,k);
324			end
325			end
326	jmc	1.27
327	mlosch	1.20	% Attributes for structure
328			if iters>0; S.iters_from_file=iters; end
329			S.attributes.global=ncgetatt(nc,'global');
330			[pstr,netcdf_fname,ext] = fileparts(fname);
331			if strcmp(netcdf_fname(end-3:end),'glob')
332			% assume it is a global file produced by gluemnc and change some
333	jmc	1.27	% attributes
334	mlosch	1.20	S.attributes.global.sNx = S.attributes.global.Nx;
335			S.attributes.global.sNy = S.attributes.global.Ny;
336			S.attributes.global.nPx = 1;
337			S.attributes.global.nSx = 1;
338			S.attributes.global.nPy = 1;
339			S.attributes.global.nSy = 1;
340			S.attributes.global.tile_number = 1;
341			S.attributes.global.nco_openmp_thread_number = 1;
342			end
343	jmc	1.27
344	mlosch	1.20	% Read variable data
345			for ivar=1:size(varlist,2)
346	jmc	1.27
347	mlosch	1.20	cvar=char(varlist{ivar});
348			varid=ncfindvarid(nc,cvar);
349			if isempty(varid)
350			disp(['No such variable ''',cvar,''' in MNC file ',fname]);
351			continue
352			end
353	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
354			% modifications make the output field smaller, let us see, if it is
355	jmc	1.27	% robust
356	mlosch	1.25	if (isfield(S,cvar) == 0); firstiter = 1; else firstiter = 0; end
357			% end code by Bruno Deremble
358	jmc	1.27
359	mlosch	1.20	[varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
360	mlosch	1.22	% does this variable contain a record (unlimited) dimension?
361			[isrecvar,recpos] = ismember(recdim,dimids);
362	jmc	1.27
363	mlosch	1.20	% Dimensions
364			clear sizVar dims
365			for k=1:length(dimids)
366			[dims{k}, sizVar(k)] = netcdf.inqDim(nc,dimids(k));
367			end
368			nDims=length(sizVar);
369	mlosch	1.22	if isrecvar
370	mlosch	1.20	if isempty(find(fii)), error('Iters not found'); end
371			it = length(dims);
372	mlosch	1.22	if length(dimids) == 1
373			% just a time or iteration variable, this will result in a vector
374			% and requires special treatment
375			icount=1;
376			tmpdata = zeros(length(fii),1);
377			for k=1:length(fii)
378			istart = fii(k)-1;
379			tmpdata(k) = netcdf.getVar(nc,varid,istart,icount,'double');
380			end
381			else
382			% from now on we assume that the record dimension is always the
383			% last dimension. This may not always be the case
384			if recpos ~= nDims
385			error(sprintf('%s\n%s\n%s%s%i%s%i', ...
386			['The current code assumes that the record ' ...
387			'dimension is the last dimension,'], ...
388			'this is not the case for variable', cvar, ...
389			': nDims = ', nDims, ...
390			', position of recDim = ', recpos))
391			end
392			istart = zeros(1,it); % indexing starts a 0
393			icount = sizVar;
394	jmc	1.27	% we always want to get only on time slice at a time
395			icount(recpos) = 1;
396	mlosch	1.22	% make your life simpler by putting the time dimension first
397			tmpdata = zeros([length(fii) sizVar(1:end-1)]);
398			for k=1:length(fii)
399			istart(recpos) = fii(k)-1; % indexing starts at 0
400			tmp = netcdf.getVar(nc,varid,istart,icount,'double');
401			tmpdata(k,:) = tmp(:);
402			end
403			% move time dimension to the end ...
404			tmpdata = shiftdim(tmpdata,1);
405			% ... and restore original shape
406			tmpdata = reshape(tmpdata,[sizVar(1:end-1) length(fii)]);
407			end
408	mlosch	1.20	else
409			it = 0;
410			tmpdata = netcdf.getVar(nc,varid,'double');
411			end
412	mlosch	1.22	%
413	mlosch	1.20	if dBug > 1,
414			fprintf([' var:',cvar,': nDims=%i ('],nDims);fprintf(' %i',sizVar);
415	jmc	1.27	fprintf('):%i,nD=%i,it=%i ;',length(size(tmpdata)),length(dims),it);
416	mlosch	1.20	end
417			[ni nj nk nm nn no np]=size(tmpdata);
418	mlosch	1.22	%
419	mlosch	1.20	[i0,j0,fn]=findTileOffset(S);
420			cdim=dims{1}; if cdim(1)~='X'; i0=0; end
421			cdim=dims{1}; if cdim(1)=='Y'; i0=j0; j0=0; end
422			if length(dims)>1;
423			cdim=dims{2}; if cdim(1)~='Y'; j0=0; end
424			else
425			j0=0;
426			end
427			if dBug > 1,
428			fprintf(' i0,ni= %i %i; j,nj= %i %i; nk=%i :',i0,ni,j0,nj,nk);
429			end
430	mlosch	1.25	% code by Bruno Deremble: if you do not want to read all tiles these
431			% modifications make the output field smaller, let us see, if it is
432	jmc	1.27	% robust
433	mlosch	1.25	if (firstiter)
434	mlosch	1.26	S.attributes.i_first.(cvar) = i0;
435			S.attributes.j_first.(cvar) = j0;
436	jmc	1.27	end
437	mlosch	1.26	i0 = i0 - S.attributes.i_first.(cvar);
438			j0 = j0 - S.attributes.j_first.(cvar);
439	mlosch	1.25	% end code by Bruno Deremble
440
441	mlosch	1.20	Sstr = '';
442			for istr = 1:max(nDims,length(dims)),
443			if istr == it, Sstr = [Sstr,'dii,'];
444			elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
445			elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
446			elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
447			elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
448			elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
449			elseif istr == 6, Sstr = [Sstr,'(1:no),'];
450			elseif istr == 7, Sstr = [Sstr,'(1:np),'];
451			else, error('Can''t handle this many dimensions!');
452			end
453			end
454			eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
455			%S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
456			if dBug > 1, fprintf(' %i',size(S.(cvar))); fprintf('\n'); end
457	mlosch	1.22	%
458	mlosch	1.20	S.attributes.(cvar)=ncgetatt(nc,cvar);
459			% replace missing or FillValues with NaN
460			if isfield(S.attributes.(cvar),'missing_value');
461			misval = S.attributes.(cvar).missing_value;
462			S.(cvar)(S.(cvar) == misval) = NaN;
463			end
464			if isfield(S.attributes.(cvar),'FillValue_');
465			misval = S.attributes.(cvar).FillValue_;
466			S.(cvar)(S.(cvar) == misval) = NaN;
467			end
468	jmc	1.27
469	mlosch	1.20	end % for ivar
470	mlosch	1.21
471			% close the file
472			netcdf.close(nc);
473	jmc	1.27
474	mlosch	1.20	if isempty(S)
475			error('Something didn''t work!!!');
476			end
477	jmc	1.27
478	mlosch	1.20	return
479
480	mlosch	1.21	function vf = ncgetvar(fname,varname)
481	mlosch	1.20	% read a netcdf variable
482	jmc	1.27
483	mlosch	1.21	nc=netcdf.open(fname,'NC_NOWRITE');
484			% find out basics about the files
485	mlosch	1.20	[ndims nvars natts dimm] = netcdf.inq(nc);
486			vf = [];
487			varid = [];
488			for k=0:nvars-1
489			if strcmp(netcdf.inqVar(nc,k),varname)
490	mlosch	1.24	varid = netcdf.inqVarID(nc,varname);
491	mlosch	1.20	end
492			end
493	jmc	1.27	if ~isempty(varid);
494	mlosch	1.20	[varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
495			% get data
496			vf = double(netcdf.getVar(nc,varid));
497			else
498			% do nothing
499			end
500	mlosch	1.21	netcdf.close(nc);
501	jmc	1.27
502	mlosch	1.20	return
503
504			function misval = ncgetmisval(nc,varid)
505
506			[varname,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
507			misval = [];
508			for k=0:natts-1
509			attn = netcdf.inqAttName(nc,varid,k);
510			if strcmp(attn,'missing_value') \| strcmp(attn,'_FillValue')
511			misval = double(netcdf.getAtt(nc,varid,attname));
512			end
513			end
514	jmc	1.27
515	mlosch	1.20	function A = ncgetatt(nc,varname)
516			% get all attributes and put them into a struct
517	jmc	1.27
518	mlosch	1.20	% 1. get global properties of file
519			[ndims nvars natts dimm] = netcdf.inq(nc);
520
521			% get variable ID and properties
522			if strcmp('global',varname)
523			% "variable" is global
524			varid = netcdf.getConstant('NC_GLOBAL');
525			else
526			% find variable ID and properties
527			varid = ncfindvarid(nc,varname);
528			if ~isempty(varid)
529			[varn,xtype,dimids,natts] = netcdf.inqVar(nc,varid);
530			else
531			warning(sprintf('variable %s not found',varname))
532			end
533			end
534
535			if natts > 1
536			for k=0:natts-1
537			attn = netcdf.inqAttName(nc,varid,k);
538			[xtype attlen]=netcdf.inqAtt(nc,varid,attn);
539			attval = netcdf.getAtt(nc,varid,attn);
540			if ~ischar(attval)
541			attval = double(attval);
542			end
543	mlosch	1.23	if strcmp(attn,'_FillValue')
544			% matlab does not allow variable names to begin with an
545			% underscore ("_"), so we have to do change the name of this
546			% obsolete attribute.
547			A.FillValue_=attval;
548			else
549			A.(char(attn))=attval;
550			end
551	mlosch	1.20	end
552			else
553			A = 'none';
554			end
555	jmc	1.27
556	mlosch	1.20	return
557
558	jmc	1.27
559	mlosch	1.20	function varid = ncfindvarid(nc,varname)
560
561			[ndims nvars natts dimm] = netcdf.inq(nc);
562			varid=[];
563			for k=0:nvars-1
564			if strcmp(netcdf.inqVar(nc,k),varname);
565	mlosch	1.24	varid = netcdf.inqVarID(nc,varname);
566	mlosch	1.20	break
567			end
568			end
569
570			return