Diagnostics/DiagUtility/rdmnc_mod2.m

function [S] = rdmnc_mod(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 mutliple fields
%
% Description:
%   This function 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.
%
% Example:
%   >> S=rdmnd('state.*','XC','YC','RC','T');
%   >> imagesc( S.XC, S.YC, S.T(:,:,1)' );

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

% 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 arguement']);
        end
    end
end

% Create list of filenames
for eachfile=file
        filepathtemp=eachfile{:};
        indecies = find(filepathtemp=='/');
        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}];
        nc=netcdf(char(eachfile),'read');
        iters = [iters,nc{'iter'}(:)];
        close(nc);
    end
    iters = unique(iters);
end

% Cycle through files for data extraction.
S.attributes=[];
for ieachfile=1:length(files)
    eachfile = [filepaths{ieachfile},files{ieachfile}];
    nc=netcdf(char(eachfile),'read');
    S=rdmnc_local(nc,varlist,iters,S);
    close(nc);
end

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

function [A] = read_att(nc);
    allatt=ncnames(att(nc));
    A='none';
    for attr=allatt;
        A.(char(attr))=nc.(char(attr))(:);
    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)

    fiter = nc{'iter'}(:);                               % File iterations present
    [fii,dii] = ismember(fiter,iters);  fii = find(fii); % File iteration index
    dii = dii(find(dii ~= 0));                           % Data interation index

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

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

        % 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
        
        dims=ncnames(dim(nc{cvar}));        % Dimensions
        
        if dims{1}=='T'
            if isempty(find(fii)), return, end
            tmpdata=nc{cvar}(find(fii),:);
            it = length(dims);
        else
            tmpdata=nc{cvar}(:);
            it = 0;
        end
        
        tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
        [ni nj nk nm nn no np]=size(tmpdata);
        [ni nj nk nm nn no np];
        
        [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
        
        Sstr = '';
        for istr = 1:7
            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;
        S.attributes.(cvar)=read_att(nc{cvar});
        end

if isempty(S)
    error('Something didn''t work!!!');
end
1	enderton	1.1	function [S] = rdmnc_mod(varargin)
2
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			%
16			% Output:
17			% S Structure with mutliple fields
18			%
19			% Description:
20			% This function rudimentary wrapper for joining and reading netcdf files
21			% produced by MITgcm. It does not give the same flexibility as opening
22			% the netcdf files directly using netcdf(). It is useful for quick
23			% loading of entire model fields which are distributed in multiple netcdf
24			% files.
25			%
26			% Example:
27			% >> S=rdmnd('state.*','XC','YC','RC','T');
28			% >> imagesc( S.XC, S.YC, S.T(:,:,1)' );
29
30			% Initializations
31			file={};
32			filepaths={};
33			files={};
34			varlist={};
35			iters=[];
36
37			% Process function arguments
38			for iarg=1:nargin;
39			arg=varargin{iarg};
40			if ischar(arg)
41			if isempty(dir(char(arg)))
42			varlist{end+1}=arg;
43			else
44			file{end+1}=arg;
45			end
46			else
47			if isempty(iters)
48			iters=arg;
49			else
50			error(['The only allowed numeric argument is iterations',...
51			' to read in as a vector for the last arguement']);
52			end
53			end
54			end
55
56			% Create list of filenames
57			for eachfile=file
58			filepathtemp=eachfile{:};
59			indecies = find(filepathtemp=='/');
60			if ~isempty(indecies)
61			filepathtemp = filepathtemp(1:indecies(end));
62			else
63			filepathtemp = '';
64			end
65			expandedEachFile=dir(char(eachfile{1}));
66			for i=1:size(expandedEachFile,1);
67			if expandedEachFile(i).isdir==0
68			files{end+1}=expandedEachFile(i).name;
69			filepaths{end+1}=filepathtemp;
70			end
71			end
72			end
73
74
75			% If iterations unspecified, open all the files and make list of all the
76			% iterations that appear, use this iterations list for data extraction.
77			if isempty(iters)
78			iters = [];
79			for ieachfile=1:length(files)
80			eachfile = [filepaths{ieachfile},files{ieachfile}];
81			nc=netcdf(char(eachfile),'read');
82			iters = [iters,nc{'iter'}(:)];
83	enderton	1.2	close(nc);
84	enderton	1.1	end
85			iters = unique(iters);
86			end
87
88			% Cycle through files for data extraction.
89			S.attributes=[];
90			for ieachfile=1:length(files)
91			eachfile = [filepaths{ieachfile},files{ieachfile}];
92			nc=netcdf(char(eachfile),'read');
93			S=rdmnc_local(nc,varlist,iters,S);
94			close(nc);
95			end
96
97			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
98			% Local functions %
99			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100
101			function [A] = read_att(nc);
102			allatt=ncnames(att(nc));
103			A='none';
104			for attr=allatt;
105			A.(char(attr))=nc.(char(attr))(:);
106			end
107
108			function [i0,j0,fn] = findTileOffset(S);
109			fn=0;
110			if isfield(S,'attributes') & isfield(S.attributes,'global')
111			G=S.attributes.global;
112			tn=G.tile_number;
113			snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
114			ntx=nsxnpx;nty=nsynpy;
115			gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
116			i0=snxgbi; j0=snygbj;
117			if isfield(S.attributes.global,'exch2_myFace')
118			fn=G.exch2_myFace;
119			end
120			else
121			i0=0;j0=0;
122			end
123			%[snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
124
125			function [S] = rdmnc_local(nc,varlist,iters,S)
126
127			fiter = nc{'iter'}(:); % File iterations present
128			[fii,dii] = ismember(fiter,iters); fii = find(fii); % File iteration index
129			dii = dii(find(dii ~= 0)); % Data interation index
130
131			% No variables specified? Default to all
132			if isempty(varlist), varlist=ncnames(var(nc)); end
133
134			% Attributes for structure
135			if iters>0; S.iters_read_from_file=iters; end
136			S.attributes.global=read_att(nc);
137
138			% Read variable data
139			for ivar=1:size(varlist,2);
140			cvar=char(varlist{ivar});
141			if isempty(nc{cvar})
142			disp(['No such variable ''',cvar,''' in MNC file ',name(nc)]);
143			continue
144			end
145
146			dims=ncnames(dim(nc{cvar})); % Dimensions
147
148			if dims{1}=='T'
149			if isempty(find(fii)), return, end
150			tmpdata=nc{cvar}(find(fii),:);
151			it = length(dims);
152			else
153			tmpdata=nc{cvar}(:);
154			it = 0;
155			end
156
157			tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
158			[ni nj nk nm nn no np]=size(tmpdata);
159			[ni nj nk nm nn no np];
160
161			[i0,j0,fn]=findTileOffset(S);
162			cdim=dims{end}; if cdim(1)~='X'; i0=0; end
163			cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
164			if length(dims)>1;
165			cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
166			else
167			j0=0;
168			end
169
170			Sstr = '';
171			for istr = 1:7
172			if istr == it, Sstr = [Sstr,'dii,'];
173			elseif istr == 1, Sstr = [Sstr,'i0+(1:ni),'];
174			elseif istr == 2, Sstr = [Sstr,'j0+(1:nj),'];
175			elseif istr == 3, Sstr = [Sstr,'(1:nk),'];
176			elseif istr == 4, Sstr = [Sstr,'(1:nm),'];
177			elseif istr == 5, Sstr = [Sstr,'(1:nn),'];
178			elseif istr == 6, Sstr = [Sstr,'(1:no),'];
179			elseif istr == 7, Sstr = [Sstr,'(1:np),'];
180			else, error('Can''t handle this many dimensions!');
181			end
182			end
183			eval(['S.(cvar)(',Sstr(1:end-1),')=tmpdata;'])
184			%S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
185			S.attributes.(cvar)=read_att(nc{cvar});
186			end
187
188			if isempty(S)
189			error('Something didn''t work!!!');
190			end