/[MITgcm]/MITgcm/utils/matlab/mnc_assembly.m
ViewVC logotype

Annotation of /MITgcm/utils/matlab/mnc_assembly.m

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


Revision 1.1 - (hide annotations) (download)
Tue Feb 1 15:55:52 2005 UTC (19 years, 3 months ago) by edhill
Branch: MAIN
 o initial checkin

1 edhill 1.1 function [nt,nf] = mnc_assembly(fpat,vars, fout,fsize)
2    
3     % Function [nt,nf] = mnc_assembly(fpat,vars, fout,fsize)
4     %
5     % INPUTS
6     % fpat string containing the file pattern
7     % vars cell array of variable names
8     %
9     % fout output file pattern (DEF: "all.%05d.nc")
10     % fsize max output file size (DEF: 2.0e+9 = +/-2GB)
11     %
12     % OUTPUTS
13     % nt number of usable tiles found
14     % nf number of output files written
15     %
16     % This function "assembles" MNC output. It finds all the per-tile
17     % NetCDF files that match the input pattern, does some basic "sanity"
18     % tests to determine whether the files have compatible sizes, and
19     % then assembles all of the requested data (all of the variables)
20     % into one or more "global" NetCDF files. The global files have
21     % the following dimension conventions:
22     %
23     % "exch 1": all values are within a global horizontal grid
24     % and indicies are (X,Y,Z,T)
25     %
26     % "exch 2": all values are within one of up to six "faces"
27     % of a global cube with indicies (Xf,Yf,F,Z,T)
28     %
29     % where "X,Y.Z,T" are global space/time indicies, "Xf,Yf" are local
30     % per-face spatial indicies, and "F" is a face index.
31    
32    
33     %===== Argument checking and defaults =====
34    
35     if nargin < 2
36     disp('Error: there must be at least 2 arguments!');
37     return
38     end
39    
40     if nargin < 3
41     fout = 'all.%05d.nc';
42     end
43     if nargin < 4
44     fsize = 2.0e+9;
45     end
46    
47    
48     %===== Find and open all the matching files =====
49    
50     nt = 0;
51     nf = 0;
52     all_ncf = struct([]);
53    
54     % Find all of the files
55     exch2_msg = 0;
56     tmax = 200;
57     frdone = 0;
58     it = 0;
59     while frdone == 0
60    
61     it = it + 1;
62     fnm = sprintf(fpat,it);
63     % disp(fnm);
64    
65     % Check that the file exists
66     fid = fopen(fnm, 'r');
67     if fid < 0
68     if it >= tmax
69     frdone = 1;
70     end
71     continue;
72     end
73    
74     % Open the NetCDF file
75     fnc = netcdf(fnm, 'nowrite');
76     if length(fnc) == 0
77     continue;
78     end
79    
80     % Check for exch1/exch2 grid
81     exch = 1;
82     exch2_myFace = fnc.exch2_myFace(:);
83     if length(exch2_myFace) ~= 0
84     exch = 2;
85     if exch2_msg == 0
86     exch2_msg = 1;
87     disp(' Grid type appears to be: "exch2"');
88     end
89     end
90    
91     n = length(all_ncf) + 1;
92     all_ncf(n).name = fnm;
93     all_ncf(n).nc = {fnc};
94     all_ncf(n).exch = exch;
95     all_ncf(n).tile_number = fnc.tile_number(1);
96     all_ncf(n).bi = fnc.bi(1);
97     all_ncf(n).bj = fnc.bj(1);
98     all_ncf(n).sNx = fnc.sNx(1);
99     all_ncf(n).sNy = fnc.sNy(1);
100     all_ncf(n).Z = fnc.Z(1);
101    
102     if exch == 2
103     all_ncf(n).exch2_myFace = exch2_myFace;
104     all_ncf(n).exch2_tbasex = fnc.exch2_tbasex(1);
105     all_ncf(n).exch2_tbasey = fnc.exch2_tbasex(1);
106     end
107    
108     clear fnc
109     end
110    
111    
112     %===== Do some basic sanity checks =====
113    
114     % check for number of files/tiles found
115     if length(all_ncf) == 0
116     disp('Error: no tiles found--no need to do any assembly!');
117     return
118     elseif length(all_ncf) == 1
119     disp('Error: one tile found--no need to do any assembly!');
120     return
121     else
122     disp(sprintf(' Found %d files matching the pattern: "%s"', ...
123     length(all_ncf), fpat ));
124     end
125    
126     % check for consistent "exch" version
127     if prod(double([all_ncf.exch] == all_ncf(1).exch)) ~= 1
128     disp('Error: not all the "exch" types of the files match.');
129     return;
130     end
131    
132     % check for consistent sNx,sNy
133     if (prod(double([all_ncf.sNx] == all_ncf(1).sNx)) ~= 1) ...
134     || (prod(double([all_ncf.sNy] == all_ncf(1).sNy)) ~= 1)
135     disp('Error: the "sNx,sNy" values for all the tiles are not');
136     disp(' uniform. Future versions of this function will be');
137     disp(' able to handle non-uniform grid sizes but this');
138     disp(' feature is not yet implemented.');
139     return;
140     end
141    
142     % check for redundant tiles and "time series" output
143     if length(all_ncf) ~= length(unique([all_ncf.tile_number]))
144     disp('Error: redundant tiles were found. Please check that');
145     disp(' the file pattern does not specify output spanning');
146     disp(' multiple model runs or even multiple time series');
147     disp(' within a single model run. For multi-time-series');
148     disp(' data sets, EACH "LEVEL" IN THE OUTPUT SERIES MUST');
149     disp(' BE ASSEMBLED SEPARATERLY.');
150     return
151     end
152    
153    
154     %===== Get the dims/vars associations =====
155    
156     mydims = struct('names', {}, 'lens', {});
157     myvars = struct([]);
158     clear tncf;
159     for ivar = 1:length(vars)
160     mydim_names = {};
161     mydim_sizes = {};
162     myatt.names = {};
163     myatt.types = {};
164     myatt.data = {};
165     myname = vars(ivar).name;
166     disp([' Looking for variable: ' myname]);
167    
168     itile = 1;
169     tncf = all_ncf(itile).nc{1};
170     ncv = tncf{myname};
171     len = length(ncv);
172     if length(ncv) == 0
173     warns = [' Warning: variable "%s" is not defined in "%s"\n' ...
174     ' so it will be ignored.'];
175     disp(sprintf(warns,myname,all_ncf(itile).name));
176     continue
177     end
178     mytype = datatype(ncv);
179     tmpdims = dim(ncv);
180     for inm = 1:length(tmpdims)
181     mydim_names{inm} = name(tmpdims{inm});
182     mydim_sizes{inm} = tmpdims{inm}(:);
183     end
184     for iat = 1:length(att(ncv))
185     aaa = att(ncv);
186     myatt.names(iat) = { name(aaa{iat}) };
187     myatt.types(iat) = { datatype(aaa{iat}) };
188     aab = aaa{iat};
189     myatt.data(iat) = { aab(:) };
190     end
191    
192     % confirm: vars have same dim names across all files
193     ierr = 0;
194     for itile = 2:length(all_ncf)
195     tncf = all_ncf(itile).nc{1};
196     ncv = tncf{myname};
197     len = length(ncv);
198     if length(ncv) == 0
199     warns = [' Warning: variable "%s" is not defined in "%s"\n' ...
200     ' so it will be ignored.'];
201     disp(sprintf(warns,myname,all_ncf(itile).name));
202     continue
203     end
204     tmpdims = dim(ncv);
205     for inm = 1:length(tmpdims)
206     if mydim_names{inm} ~= name(tmpdims{inm})
207     warns = ...
208     [' Warning: variable "%s" is not CONSISTENTLY defined.\n' ...
209     ' It has different dimensions in different files so\n' ...
210     ' so it will be ignored.'];
211     disp(sprintf(warns,myname));
212     ierr = 1;
213     break
214     end
215     mydim_sizes{inm} = max([ tmpdims{inm}(:) mydim_sizes{inm} ]);
216     end
217    
218     end
219    
220     if ierr == 0
221     % check: does the variable have a "horizontal" component
222     has_horiz = 0;
223     horiz_names = { 'X' 'Y' 'Xp1' 'Yp1' };
224     for id = 1:length(mydim_names)
225     if length([intersect(horiz_names,mydim_names{id})]) > 0
226     has_horiz = 1;
227     end
228     end
229     % disp([ ' ' myname ' ' sprintf('%d',has_horiz) ]);
230    
231     imy = length(myvars) + 1;
232     myvars(imy).name = myname;
233     myvars(imy).type = mytype;
234     myvars(imy).dim_names = mydim_names;
235     myvars(imy).dim_sizes = mydim_sizes;
236     myvars(imy).atts = myatt;
237     myvars(imy).has_horiz = has_horiz;
238    
239     addl = setdiff(mydim_names,[mydims.names]);
240     for iaddl = 1:length(addl)
241     np1 = length(mydims) + 1;
242     mydims(np1).names = addl(iaddl);
243     mydims(np1).lens = mydim_sizes(find(strcmp(addl(iaddl),mydim_names)));
244     end
245    
246     end
247     end
248    
249     % For exch == 2, we need to add a "face" dimension
250     if all_ncf(1).exch == 2
251     np1 = length(mydims) + 1;
252     mydims(np1).names = { 'iface' };
253     mydims(np1).lens = { length(unique([all_ncf.exch2_myFace])) };
254     end
255    
256     % myvars.name
257     % myvars.dim_names
258     % myvars.dim_sizes
259     % myvars(2).dim_names
260     % myvars(2).dim_names(4)
261    
262     % mydims
263     % length(mydims)
264     % [ mydims.names ]
265     % [ mydims.lens ]
266    
267    
268     %===== Assemble! =====
269    
270    
271     if all_ncf(1).exch == 1
272    
273     % exch "1":
274    
275     bi_max = max([all_ncf.bi]);
276     bj_max = max([all_ncf.bj]);
277     Xmax = bi_max * all_ncf(1).sNx;
278     Ymax = bj_max * all_ncf(1).sNy;
279    
280     horzdim = struct('names',{},'lens',{});
281     horzdim(1).names = { 'X' }; horzdim(1).lens = { Xmax };
282     horzdim(2).names = {'Xp1'}; horzdim(2).lens = { Xmax + 1 };
283     horzdim(3).names = { 'Y' }; horzdim(3).lens = { Ymax };
284     horzdim(4).names = {'Yp1'}; horzdim(4).lens = { Ymax + 1 };
285     horzdim(5).names = { 'T' }; horzdim(5).lens = { 0 };
286    
287     iseq = 0;
288     foutnm = sprintf(fout, iseq);
289     fonc = netcdf(foutnm,'clobber'); % Should append-or-create!
290    
291     for idim = 1:length(mydims)
292     dname = mydims(idim).names{1};
293     ind = find(strcmp(dname,[horzdim.names]));
294     if length(ind) ~= 0
295     dlen = horzdim(ind).lens{1};
296     else
297     dlen = mydims(idim).lens{1};
298     end
299     comm = sprintf('fonc(''%s'') = %d;',dname,dlen);
300     eval(comm);
301     end
302    
303     for ivar = 1:length(myvars)
304     comm = sprintf('fonc{''%s''} = nc%s( ',myvars(ivar).name,myvars(ivar).type);
305     id = 1;
306     comm = [ comm sprintf('''%s''',myvars(ivar).dim_names{id}) ];
307     for id = 2:length(myvars(ivar).dim_names)
308     comm = [ comm sprintf(',''%s''',myvars(ivar).dim_names{id}) ];
309     end
310     comm = [ comm ' );' ];
311     eval(comm);
312     for iat = 1:length(myvars(ivar).atts.names)
313     comm = sprintf( ...
314     'fonc{''%s''}.%s = nc%s( myvars(ivar).atts.data{iat} );', ...
315     myvars(ivar).name, ...
316     myvars(ivar).atts.names{iat}, ...
317     myvars(ivar).atts.types{iat} );
318     eval(comm);
319     end
320     end
321    
322     % for itime = 1:Tmax
323    
324     % Here is where we need to check the output file size and start
325     % another file in the sequence, if necessary.
326    
327     for ivar = 1:length(myvars)
328     disp(sprintf(' Copying variable: %s',myvars(ivar).name))
329     for itile = 1:length(all_ncf)
330    
331     if (myvars(ivar).has_horiz == 1) || (itile == 1)
332    
333     clear nct;
334     nct = all_ncf(itile).nc{1};
335     ox_off = (all_ncf(itile).bi - 1)*all_ncf(itile).sNx;
336     oy_off = (all_ncf(itile).bj - 1)*all_ncf(itile).sNy;
337     diml_in = '';
338     diml_out = '';
339     for jj = 1:length(myvars(ivar).dim_names)
340     doff = 1;
341     if jj > 1
342     diml_in = sprintf('%s,',diml_in);
343     diml_out = sprintf('%s,',diml_out);
344     end
345     dlen = myvars(ivar).dim_sizes{jj};
346     diml_in = sprintf('%s%s',diml_in, ':');
347     fchar = myvars(ivar).dim_names{jj}(1);
348     % disp([' fchar = ' fchar ' ' myvars(ivar).dim_names{jj}]);
349     if strcmp(myvars(ivar).dim_names{jj}(1),'X') == 1
350     doff = ox_off + doff;
351     dlen = ox_off + dlen;
352     end
353     if strcmp(myvars(ivar).dim_names{jj}(1),'Y') == 1
354     doff = oy_off + doff;
355     dlen = oy_off + dlen;
356     end
357     diml_out = sprintf('%s%d%s%d',diml_out,doff,':',dlen);
358     end
359    
360     comm = sprintf( ...
361     'fonc{''%s''}(%s) = nct{''%s''}(%s);', ...
362     myvars(ivar).name, diml_out, myvars(ivar).name, diml_in );
363     % disp([ ' comm: ' comm ]);
364     eval(comm);
365    
366     end
367    
368     end
369     end
370     % end
371    
372     fonc = close(fonc);
373    
374     elseif all_ncf(1).exch == 2
375    
376     % exch "2":
377     Xmax = 0;
378     Ymax = 0;
379     for ii = 1:length(all_ncf)
380     Xmax = max(Xmax, (all_ncf(ii).exch2_tbasex + all_ncf(ii).sNx));
381     Ymax = max(Ymax, (all_ncf(ii).exch2_tbasey + all_ncf(ii).sNy));
382     end
383    
384     horzdim = struct('names',{},'lens',{});
385     horzdim(1).names = { 'X' }; horzdim(1).lens = { Xmax };
386     horzdim(2).names = {'Xp1'}; horzdim(2).lens = { Xmax + 1 };
387     horzdim(3).names = { 'Y' }; horzdim(3).lens = { Ymax };
388     horzdim(4).names = {'Yp1'}; horzdim(4).lens = { Ymax + 1 };
389     horzdim(5).names = { 'T' }; horzdim(5).lens = { 0 };
390    
391     iseq = 0;
392     foutnm = sprintf(fout, iseq);
393     fonc = netcdf(foutnm,'clobber'); % Should append-or-create!
394    
395     for idim = 1:length(mydims)
396     dname = mydims(idim).names{1};
397     ind = find(strcmp(dname,[horzdim.names]));
398     if length(ind) ~= 0
399     dlen = horzdim(ind).lens{1};
400     else
401     dlen = mydims(idim).lens{1};
402     end
403     comm = sprintf('fonc(''%s'') = %d;',dname,dlen);
404     eval(comm);
405     end
406    
407     for ivar = 1:length(myvars)
408     comm = sprintf('fonc{''%s''} = nc%s( ',myvars(ivar).name,myvars(ivar).type);
409     id = 1;
410     comm = [ comm sprintf('''%s''',myvars(ivar).dim_names{id}) ];
411     for id = 2:length(myvars(ivar).dim_names)
412     dname = myvars(ivar).dim_names{id};
413     if (dname(1) == 'Y') && (myvars(ivar).has_horiz == 1)
414     comm = [ comm sprintf(',''%s''','iface') ];
415     end
416     comm = [ comm sprintf(',''%s''',dname) ];
417     end
418     comm = [ comm ' );' ];
419     eval(comm);
420     for iat = 1:length(myvars(ivar).atts.names)
421     comm = sprintf( ...
422     'fonc{''%s''}.%s = nc%s( myvars(ivar).atts.data{iat} );', ...
423     myvars(ivar).name, ...
424     myvars(ivar).atts.names{iat}, ...
425     myvars(ivar).atts.types{iat} );
426     eval(comm);
427     end
428     end
429    
430     % Here is where we need to check the output file size and start
431     % another file in the sequence, if necessary.
432    
433     for ivar = 1:length(myvars)
434     disp(sprintf(' Copying variable: %s',myvars(ivar).name))
435     for itile = 1:length(all_ncf)
436    
437     if (myvars(ivar).has_horiz == 1) || (itile == 1)
438    
439     clear nct;
440     nct = all_ncf(itile).nc{1};
441     ox_off = all_ncf(itile).exch2_tbasex;
442     oy_off = all_ncf(itile).exch2_tbasey;
443     diml_tin = '';
444     diml_res = '';
445     diml_in = '';
446     diml_out = '';
447     if length(myvars(ivar).dim_names) < 2
448     comm = sprintf( ...
449     'fonc{''%s''}(%s%d) = nct{''%s''}(:);', ...
450     myvars(ivar).name, '1:', myvars(ivar).dim_sizes{1}, ...
451     myvars(ivar).name );
452     % disp([ ' ' comm ]);
453     eval(comm);
454     else
455     for jj = 1:length(myvars(ivar).dim_names)
456     doff = 1;
457     if jj > 1
458     diml_tin = sprintf('%s,',diml_tin);
459     diml_res = sprintf('%s,',diml_res);
460     diml_in = sprintf('%s,',diml_in);
461     diml_out = sprintf('%s,',diml_out);
462     end
463     dnam = myvars(ivar).dim_names{jj};
464     dlen = myvars(ivar).dim_sizes{jj};
465     dlenr = dlen;
466     fchar = myvars(ivar).dim_names{jj}(1);
467     % disp([' fchar = ' fchar ' ' myvars(ivar).dim_names{jj}]);
468     if strcmp(dnam(1),'X') == 1
469     doff = ox_off + doff;
470     dlen = ox_off + dlen;
471     end
472     if strcmp(dnam(1),'Y') == 1
473     diml_res = sprintf('%s%s',diml_res, '[],');
474     diml_in = sprintf('%s%s',diml_in, ':,');
475     diml_out = sprintf('%s%d%s',diml_out,all_ncf(itile).exch2_myFace,',');
476     doff = oy_off + doff;
477     dlen = oy_off + dlen;
478     end
479     diml_tin = sprintf('%s%s',diml_tin, ':');
480     diml_res = sprintf('%s%d',diml_res, dlenr);
481     diml_in = sprintf('%s%s',diml_in, ':');
482     diml_out = sprintf('%s%d%s%d',diml_out,doff,':',dlen);
483     end
484    
485     comm = sprintf( ...
486     'tmp = reshape(nct{''%s''}(%s), %s); fonc{''%s''}(%s) = tmp(%s);', ...
487     myvars(ivar).name, diml_tin, diml_res, myvars(ivar).name, ...
488     diml_out, diml_in );
489     % disp([ ' ' comm ]);
490     eval(comm);
491     end
492    
493     end
494    
495     end
496     end
497     % end
498    
499     fonc = close(fonc);
500    
501     end
502    

  ViewVC Help
Powered by ViewVC 1.1.22