utils/matlab/rdmnc.m

function [S] = rdmnc(varargin)
% S=RDMNC(FILE1,FILE2,...)
% S=RDMNC(FILE1,...,ITER)
% S=RDMNC(FILE1,...,'VAR1','VAR2',...)
% S=RDMNC(FILE1,...,'VAR1','VAR2',...,ITER)
%
%  FILE1 ... is 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 ... are model variable names as written in the MNC/netcdf file
%  ITER  is a vector of time levels (consecutive indices referring to snap-shot
%        in the MNC/netcdf file i.e. 1,2,3,... and not model time)
%  S     is a structure with mutliple fields
%
% rdmnc() 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(). It is useful for quick loading of
% entire model fields which are distributed in multiple netcdf files.
% rdmnc() also reads non-MITgcm generated netcdf files.
%
% e.g. To plot the surface temperature in last time-level written to file
% >> S=rdmnd('state.*','XC','YC','RC','T',Inf);
% >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
%
% $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.1 2004/03/08 16:45:56 adcroft Exp $

%Defaults
global verbose
file={};
files={};
verbose=0;
varlist={};
timelevels=[];

% Process function arguments
for iarg=1:nargin;
 arg=varargin{iarg};
 if ischar(arg)
   switch char(arg)
     case 'verbose',
       verbose=1;
     otherwise,
       if isempty( dir(char(arg)) )
        varlist{end+1}={arg};
       else
        file{end+1}={arg};
       end
   end
 else
   if isempty(timelevels)
     timelevels=arg;
   else
     error('Specify the time levels in a vector and not as multiple numeric arguments')
   end
 end
end

if verbose; disp('Verbose mode enabled'); end

% Create list of filenames
for eachfile=file
 expandedEachFile=dir(char(eachfile{1}));
 for i=1:size(expandedEachFile,1);
  if expandedEachFile(i).isdir==0; files{end+1}=expandedEachFile(i).name; end
 end
end

% Open each file
S.attributes=[];
for eachfile={files{1:end}}
 nc=netcdf(char(eachfile),'read');
 if ismncfile(nc)
  S=rdmnc_local(nc,varlist,timelevels,S);
 else
  S=rdnetcdf_local(nc,varlist,S);
 end
end

% -----------------------------------------------------------------------------
function [result] = ismncfile(nc);
result=~isempty(nc.MITgcm_mnc_ver);
%MLresult=~isempty(nc.('MITgcm_mnc_ver'));
% -----------------------------------------------------------------------------
function [A] = read_att(nc);
global verbose
allatt=ncnames(att(nc)); if verbose; allatt, end
A='none';
for attr=allatt;
 tmp = char(attr);
 eval(['A.' tmp '= nc.' tmp '(:);'])
%ML A.(char(attr))=nc.(char(attr))(:);
end
% -----------------------------------------------------------------------------
function [i0,j0,fn] = findTileOffset(S);
global verbose
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
% -----------------------------------------------------------------------------
function [S] = rdnetcdf_local(nc,varlist,S)
% Read all attributes and variable data from a netcdf file
% with special operations for MNC style data
global verbose

% No variables specified? Default to all
if isempty(varlist)
 varlist=ncnames(var(nc)); if verbose; varlist, end
end

% Attributes for structure
S.attributes=read_att(nc);

% Read variable data
for ivar=1:size(varlist,2);
 cvar=char(varlist{ivar});
 tmpdata=nc{cvar}(:);
 if isempty(tmpdata)
  disp(['No such variable ''' cvar ''' in netcdf file' name(nc)])
 else
  tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
  eval(['S.' cvar '=tmpdata;'])
  eval(['S.attributes.' cvar '=read_att(nc{' cvar '});'])
%ML  S.(cvar)=tmpdata;
%ML  S.attributes.(cvar)=read_att(nc{cvar});
 end
end
% -----------------------------------------------------------------------------
function [S] = rdmnc_local(nc,varlist,timelevels,S)
% Read all attributes and variable data from a netcdf file
% with special operations for MNC style data
global verbose

nt=size( nc('T'), 1); if verbose; nt, end

% No variables specified? Default to all
if isempty(varlist)
 varlist=ncnames(var(nc)); if verbose; varlist, end
end

% No iterations specified? Default to all
if isempty(timelevels) | isnan(timelevels)
 timelevels=1:nt;
elseif timelevels == Inf
 timelevels=nt;
end

% Sanity checks
if max( find(timelevels > nt) )
 error('Requested time level is beyond time dimension in netcdf file')
end

% Attributes for structure
if timelevels>0; S.timelevels_read_from_file=timelevels; end
S.attributes.global=read_att(nc);

% Read variable data
for ivar=1:size(varlist,2);
 cvar=char(varlist{ivar}); if verbose; cvar, end
 if isempty(nc{cvar})
  disp(['No such variable ''' cvar ''' in MNC file ' name(nc)])
  continue
 end
 dims=ncnames(dim(nc{cvar}));
 if dims{1}=='T'
  if verbose; disp(['Reading variable ''' cvar ''' with record dimension']); end
  tmpdata=nc{cvar}(timelevels,:);
 else
  if verbose; disp(['Reading variable ''' cvar '''']); end
  tmpdata=nc{cvar}(:);
 end
 if isempty(tmpdata)
  error(['No such variable ''' cvar ''' in MNC file ' name(nc)])
 else
  tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
  [ni nj nk nm nn no np]=size(tmpdata);
  if np~=1; error('Wow! This is a very high dimension variable...'); end
  [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
  eval(['S.' cvar ...
        '(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;'])
  eval(['S.attributes.' cvar ' =read_att(nc{''' cvar '''});'])
%ML  S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
%ML  S.attributes.(cvar)=read_att(nc{cvar});
 end
end

if isempty(S)
 error('Something didn''t work!!!')
end
% -----------------------------------------------------------------------------
1	adcroft	1.1	function [S] = rdmnc(varargin)
2			% S=RDMNC(FILE1,FILE2,...)
3			% S=RDMNC(FILE1,...,ITER)
4			% S=RDMNC(FILE1,...,'VAR1','VAR2',...)
5			% S=RDMNC(FILE1,...,'VAR1','VAR2',...,ITER)
6			%
7			% FILE1 ... is either a single file name (e.g. 'state.nc') or a wild-card
8			% strings expanding to a group of file names (e.g. 'state.*.nc').
9			% There are no assumptions about suffices (e.g. 'state.*' works)
10			% VAR1 ... are model variable names as written in the MNC/netcdf file
11			% ITER is a vector of time levels (consecutive indices referring to snap-shot
12			% in the MNC/netcdf file i.e. 1,2,3,... and not model time)
13			% S is a structure with mutliple fields
14			%
15			% rdmnc() is a rudimentary wrapper for joining and reading netcdf files
16			% produced by MITgcm. It does not give the same flexibility as opening the
17			% netcdf files directly using netcdf(). It is useful for quick loading of
18			% entire model fields which are distributed in multiple netcdf files.
19			% rdmnc() also reads non-MITgcm generated netcdf files.
20			%
21			% e.g. To plot the surface temperature in last time-level written to file
22			% >> S=rdmnd('state.*','XC','YC','RC','T',Inf);
23			% >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
24			%
25	mlosch	1.2	% $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.1 2004/03/08 16:45:56 adcroft Exp $
26	adcroft	1.1
27			%Defaults
28			global verbose
29			file={};
30			files={};
31			verbose=0;
32			varlist={};
33			timelevels=[];
34
35			% Process function arguments
36			for iarg=1:nargin;
37			arg=varargin{iarg};
38			if ischar(arg)
39			switch char(arg)
40			case 'verbose',
41			verbose=1;
42			otherwise,
43			if isempty( dir(char(arg)) )
44			varlist{end+1}={arg};
45			else
46			file{end+1}={arg};
47			end
48			end
49			else
50			if isempty(timelevels)
51			timelevels=arg;
52			else
53			error('Specify the time levels in a vector and not as multiple numeric arguments')
54			end
55			end
56			end
57
58			if verbose; disp('Verbose mode enabled'); end
59
60			% Create list of filenames
61			for eachfile=file
62			expandedEachFile=dir(char(eachfile{1}));
63			for i=1:size(expandedEachFile,1);
64			if expandedEachFile(i).isdir==0; files{end+1}=expandedEachFile(i).name; end
65			end
66			end
67
68			% Open each file
69			S.attributes=[];
70			for eachfile={files{1:end}}
71			nc=netcdf(char(eachfile),'read');
72			if ismncfile(nc)
73			S=rdmnc_local(nc,varlist,timelevels,S);
74			else
75			S=rdnetcdf_local(nc,varlist,S);
76			end
77			end
78
79			% -----------------------------------------------------------------------------
80			function [result] = ismncfile(nc);
81	mlosch	1.2	result=~isempty(nc.MITgcm_mnc_ver);
82			%MLresult=~isempty(nc.('MITgcm_mnc_ver'));
83	adcroft	1.1	% -----------------------------------------------------------------------------
84			function [A] = read_att(nc);
85			global verbose
86			allatt=ncnames(att(nc)); if verbose; allatt, end
87			A='none';
88			for attr=allatt;
89	mlosch	1.2	tmp = char(attr);
90			eval(['A.' tmp '= nc.' tmp '(:);'])
91			%ML A.(char(attr))=nc.(char(attr))(:);
92	adcroft	1.1	end
93			% -----------------------------------------------------------------------------
94			function [i0,j0,fn] = findTileOffset(S);
95			global verbose
96			fn=0;
97			if isfield(S,'attributes') & isfield(S.attributes,'global')
98			G=S.attributes.global;
99			tn=G.tile_number;
100			snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
101			ntx=nsxnpx;nty=nsynpy;
102			gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
103			i0=snxgbi; j0=snygbj;
104			if isfield(S.attributes.global,'exch2_myFace')
105			fn=G.exch2_myFace;
106			end
107			else
108			i0=0;j0=0;
109			end
110			% -----------------------------------------------------------------------------
111			function [S] = rdnetcdf_local(nc,varlist,S)
112			% Read all attributes and variable data from a netcdf file
113			% with special operations for MNC style data
114			global verbose
115
116			% No variables specified? Default to all
117			if isempty(varlist)
118			varlist=ncnames(var(nc)); if verbose; varlist, end
119			end
120
121			% Attributes for structure
122			S.attributes=read_att(nc);
123
124			% Read variable data
125			for ivar=1:size(varlist,2);
126			cvar=char(varlist{ivar});
127			tmpdata=nc{cvar}(:);
128			if isempty(tmpdata)
129			disp(['No such variable ''' cvar ''' in netcdf file' name(nc)])
130			else
131			tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
132	mlosch	1.2	eval(['S.' cvar '=tmpdata;'])
133			eval(['S.attributes.' cvar '=read_att(nc{' cvar '});'])
134			%ML S.(cvar)=tmpdata;
135			%ML S.attributes.(cvar)=read_att(nc{cvar});
136	adcroft	1.1	end
137			end
138			% -----------------------------------------------------------------------------
139			function [S] = rdmnc_local(nc,varlist,timelevels,S)
140			% Read all attributes and variable data from a netcdf file
141			% with special operations for MNC style data
142			global verbose
143
144			nt=size( nc('T'), 1); if verbose; nt, end
145
146			% No variables specified? Default to all
147			if isempty(varlist)
148			varlist=ncnames(var(nc)); if verbose; varlist, end
149			end
150
151			% No iterations specified? Default to all
152			if isempty(timelevels) \| isnan(timelevels)
153			timelevels=1:nt;
154			elseif timelevels == Inf
155			timelevels=nt;
156			end
157
158			% Sanity checks
159			if max( find(timelevels > nt) )
160			error('Requested time level is beyond time dimension in netcdf file')
161			end
162
163			% Attributes for structure
164			if timelevels>0; S.timelevels_read_from_file=timelevels; end
165			S.attributes.global=read_att(nc);
166
167			% Read variable data
168			for ivar=1:size(varlist,2);
169			cvar=char(varlist{ivar}); if verbose; cvar, end
170			if isempty(nc{cvar})
171			disp(['No such variable ''' cvar ''' in MNC file ' name(nc)])
172			continue
173			end
174			dims=ncnames(dim(nc{cvar}));
175			if dims{1}=='T'
176			if verbose; disp(['Reading variable ''' cvar ''' with record dimension']); end
177			tmpdata=nc{cvar}(timelevels,:);
178			else
179			if verbose; disp(['Reading variable ''' cvar '''']); end
180			tmpdata=nc{cvar}(:);
181			end
182			if isempty(tmpdata)
183			error(['No such variable ''' cvar ''' in MNC file ' name(nc)])
184			else
185			tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
186			[ni nj nk nm nn no np]=size(tmpdata);
187			if np~=1; error('Wow! This is a very high dimension variable...'); end
188			[i0,j0,fn]=findTileOffset(S);
189			cdim=dims{end}; if cdim(1)~='X'; i0=0; end
190			cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
191			if length(dims)>1;
192			cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
193			else
194			j0=0;
195			end
196	mlosch	1.2	eval(['S.' cvar ...
197			'(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;'])
198			eval(['S.attributes.' cvar ' =read_att(nc{''' cvar '''});'])
199			%ML S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
200			%ML S.attributes.(cvar)=read_att(nc{cvar});
201	adcroft	1.1	end
202			end
203
204			if isempty(S)
205			error('Something didn''t work!!!')
206			end
207			% -----------------------------------------------------------------------------