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

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

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


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

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