/[MITgcm]/MITgcm_contrib/enderton/Diagnostics/DiagUtility/rdmnc_mod.m
ViewVC logotype

Annotation of /MITgcm_contrib/enderton/Diagnostics/DiagUtility/rdmnc_mod.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Mon Jan 31 15:43:29 2005 UTC (20 years, 5 months ago) by enderton
Branch: MAIN
CVS Tags: HEAD
 o Initial check in.

1 enderton 1.1 function [S] = rdmnc_mod(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     % $Header: /u/gcmpack/MITgcm/utils/matlab/rdmnc.m,v 1.1 2004/03/08 16:45:56 adcroft Exp $
26    
27     %Defaults
28     global verbose
29     file={};
30     filepaths={};
31     files={};
32     verbose=0;
33     varlist={};
34     timelevels=[];
35    
36     % Process function arguments
37     for iarg=1:nargin;
38     arg=varargin{iarg};
39     if ischar(arg)
40     switch char(arg)
41     case 'verbose',
42     verbose=1;
43     otherwise,
44     if isempty( dir(char(arg)) )
45     varlist{end+1}={arg};
46     else
47     file{end+1}={arg};
48     end
49     end
50     else
51     if isempty(timelevels)
52     timelevels=arg;
53     else
54     error('Specify the time levels in a vector and not as multiple numeric arguments')
55     end
56     end
57     end
58    
59     if verbose; disp('Verbose mode enabled'); end
60    
61     % Create list of filenames
62     for eachfile=file
63     filepathtemp=eachfile{:}{:};
64     indecies = find(filepathtemp=='/');
65     if ~isempty(indecies)
66     filepathtemp = filepathtemp(1:indecies(end));
67     else
68     filepathtemp = '';
69     end
70     expandedEachFile=dir(char(eachfile{1}));
71     for i=1:size(expandedEachFile,1);
72     if expandedEachFile(i).isdir==0
73     files{end+1}=expandedEachFile(i).name;
74     filepaths{end+1}=filepathtemp;
75     end
76     end
77     end
78    
79     % Open each file
80     S.attributes=[];
81     for ieachfile=1:length(files)
82     eachfile = [filepaths{ieachfile},files{ieachfile}];
83     nc=netcdf(char(eachfile),'read');
84     if ismncfile(nc)
85     S=rdmnc_local(nc,varlist,timelevels,S);
86     else
87     S=rdnetcdf_local(nc,varlist,S);
88     end
89     close(nc);
90     end
91    
92     % -----------------------------------------------------------------------------
93     function [result] = ismncfile(nc);
94     result=~isempty(nc.('MITgcm_mnc_ver'));
95     % -----------------------------------------------------------------------------
96     function [A] = read_att(nc);
97     global verbose
98     allatt=ncnames(att(nc)); if verbose; allatt, end
99     A='none';
100     for attr=allatt;
101     A.(char(attr))=nc.(char(attr))(:);
102     end
103     % -----------------------------------------------------------------------------
104     function [i0,j0,fn] = findTileOffset(S);
105     global verbose
106     fn=0;
107     if isfield(S,'attributes') & isfield(S.attributes,'global')
108     G=S.attributes.global;
109     tn=G.tile_number;
110     snx=G.sNx; sny=G.sNy; nsx=G.nSx; nsy=G.nSy; npx=G.nPx; npy=G.nPy;
111     ntx=nsx*npx;nty=nsy*npy;
112     gbi=mod(tn-1,ntx); gbj=(tn-gbi-1)/ntx;
113     i0=snx*gbi; j0=sny*gbj;
114     if isfield(S.attributes.global,'exch2_myFace')
115     fn=G.exch2_myFace;
116     end
117     else
118     i0=0;j0=0;
119     end
120     [snx,sny,nsx,nsy,npx,npy,ntx,nty,i0,j0,fn];
121     % -----------------------------------------------------------------------------
122     function [S] = rdnetcdf_local(nc,varlist,S)
123     % Read all attributes and variable data from a netcdf file
124     % with special operations for MNC style data
125     global verbose
126    
127     % No variables specified? Default to all
128     if isempty(varlist)
129     varlist=ncnames(var(nc)); if verbose; varlist, end
130     end
131    
132     % Attributes for structure
133     S.attributes=read_att(nc);
134    
135     % Read variable data
136     for ivar=1:size(varlist,2);
137     cvar=char(varlist{ivar});
138     tmpdata=nc{cvar}(:);
139     if isempty(tmpdata)
140     disp(['No such variable ''' cvar ''' in netcdf file' name(nc)])
141     else
142     tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
143     S.(cvar)=tmpdata;
144     S.attributes.(cvar)=read_att(nc{cvar});
145     end
146     end
147     % -----------------------------------------------------------------------------
148     function [S] = rdmnc_local(nc,varlist,timelevels,S)
149     % Read all attributes and variable data from a netcdf file
150     % with special operations for MNC style data
151     global verbose
152    
153     nt=size( nc('T'), 1); if verbose; nt, end
154    
155     % No variables specified? Default to all
156     if isempty(varlist)
157     varlist=ncnames(var(nc)); if verbose; varlist, end
158     end
159    
160     % No iterations specified? Default to all
161     if isempty(timelevels) | isnan(timelevels)
162     timelevels=1:nt;
163     elseif timelevels == Inf
164     timelevels=nt;
165     end
166    
167     % Sanity checks
168     if max( find(timelevels > nt) )
169     error('Requested time level is beyond time dimension in netcdf file')
170     end
171    
172     % Attributes for structure
173     if timelevels>0; S.timelevels_read_from_file=timelevels; end
174     S.attributes.global=read_att(nc);
175    
176     % Read variable data
177     for ivar=1:size(varlist,2);
178     cvar=char(varlist{ivar}); if verbose; cvar, end
179     if isempty(nc{cvar})
180     disp(['No such variable ''' cvar ''' in MNC file ' name(nc)])
181     continue
182     end
183     dims=ncnames(dim(nc{cvar}));
184     if dims{1}=='T'
185     if verbose; disp(['Reading variable ''' cvar ''' with record dimension']); end
186     tmpdata=nc{cvar}(timelevels,:);
187     else
188     if verbose; disp(['Reading variable ''' cvar '''']); end
189     tmpdata=nc{cvar}(:);
190     end
191     if isempty(tmpdata)
192     error(['No such variable ''' cvar ''' in MNC file ' name(nc)])
193     else
194     tmpdata=squeeze(permute(tmpdata,[9:-1:1]));
195     [ni nj nk nm nn no np]=size(tmpdata);
196     [ni nj nk nm nn no np];
197     if np~=1; error('Wow! This is a very high dimension variable...'); end
198     [i0,j0,fn]=findTileOffset(S);
199     cdim=dims{end}; if cdim(1)~='X'; i0=0; end
200     cdim=dims{end}; if cdim(1)=='Y'; i0=j0; j0=0; end
201     if length(dims)>1;
202     cdim=dims{end-1}; if cdim(1)~='Y'; j0=0; end
203     else
204     j0=0;
205     end
206     S.(cvar)(i0+(1:ni),j0+(1:nj),(1:nk),(1:nm),(1:nn),(1:no),(1:np))=tmpdata;
207     S.attributes.(cvar)=read_att(nc{cvar});
208     end
209     end
210    
211     if isempty(S)
212     error('Something didn''t work!!!')
213     end
214     % -----------------------------------------------------------------------------

  ViewVC Help
Powered by ViewVC 1.1.22