/[MITgcm]/MITgcm_contrib/high_res_cube/eddy_flux/rdmds.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/eddy_flux/rdmds.m

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


Revision 1.1 - (hide annotations) (download)
Tue May 18 01:53:30 2004 UTC (21 years, 2 months ago) by edhill
Branch: MAIN
 o initial check-in

1 edhill 1.1 function [AA,iters] = rdmds(fnamearg,varargin)
2     % RDMDS Read MITgcmUV meta/data files
3     %
4     % A = RDMDS(FNAME)
5     % A = RDMDS(FNAME,ITER)
6     % A = RDMDS(FNAME,[ITER1 ITER2 ...])
7     % A = RDMDS(FNAME,NaN)
8     % A = RDMDS(FNAME,Inf)
9     % [A,ITS] = RDMDS(FNAME,[...])
10     % A = RDMDS(FNAME,[...],'rec',RECNUM)
11     %
12     % A = RDMDS(FNAME) reads data described by meta/data file format.
13     % FNAME is a string containing the "head" of the file names.
14     %
15     % eg. To load the meta-data files
16     % T.0000002880.000.000.meta, T.0000002880.000.000.data
17     % T.0000002880.001.000.meta, T.0000002880.001.000.data
18     % T.0000002880.002.000.meta, T.0000002880.002.000.data
19     % T.0000002880.003.000.meta, T.0000002880.003.000.data
20     % use
21     % >> A=rdmds('T.0000002880');
22     % >> size(A)
23     % ans =
24     % 64 32 5
25     % eg. To load a multiple record file
26     % >> A=rdmds('pickup.0000002880');
27     % >> size(A)
28     % ans =
29     % 64 32 5 61
30     %
31     %
32     % A = RDMDS(FNAME,ITER) reads data described by meta/data file format.
33     % FNAME is a string containing the "head" of the file name excluding the
34     % 10-digit iterartion number.
35     % ITER is a vector of positive integers that will expand to the 10-digit
36     % number in the file name.
37     % If ITER=NaN, all iterations will be read.
38     % If ITER=Inf, the last (highest) iteration will be read.
39     %
40     % eg. To repeat above operation
41     % >> A=rdmds('T',2880);
42     % eg. To read multiple time steps
43     % >> A=rdmds('T',[0 1440 2880]);
44     % eg. To read all time steps
45     % >> [A,ITS]=rdmds('T',NaN);
46     % eg. To read the last time step
47     % >> [A,IT]=rdmds('T',Inf);
48     % Note: this form can not read files with no iteration count in file name.
49     %
50     %
51     % A = RDMDS(FNAME,[...],'rec',RECNUM) reads individual records from
52     % multiple record files.
53     %
54     % eg. To read a single record from a multi-record file
55     % >> [A,IT]=rdmds('pickup.ckptA',11);
56     % eg. To read several records from a multi-record file
57     % >> [A,IT]=rdmds('pickup',Inf,'rec',[1:5 8 12]);
58     %
59     %
60     % A = RDMDS(FNAME,ITER,MACHINEFORMAT) allows the machine format to be
61     % A = RDMDS(FNAME,MACHINEFORMAT)
62     % specified which MACHINEFORMAT is on of the following strings:
63     % 'n' 'l' 'b' 'd' 'g' 'c' 'a' 's' - see FOPEN for more details
64     %
65     %
66     % $Header: /u/gcmpack/MITgcm/verification/global_ocean.cubed32x32x30/matlab/rdmds.m,v 1.1 2002/12/13 17:57:25 cheisey Exp $
67    
68     AA=[];
69     iters=[];
70    
71     % Default options
72     ieee='b';
73     fname=fnamearg;
74     userecords=0;
75     recnum=[];
76    
77     % Check optional arguments
78     for ind=1:size(varargin,2);
79     arg=varargin{ind};
80     if ischar(arg)
81     if strcmp(arg,'n') | strcmp(arg,'native')
82     ieee='n';
83     elseif strcmp(arg,'l') | strcmp(arg,'ieee-le')
84     ieee='l';
85     elseif strcmp(arg,'b') | strcmp(arg,'ieee-be')
86     ieee='b';
87     elseif strcmp(arg,'c') | strcmp(arg,'cray')
88     ieee='c';
89     elseif strcmp(arg,'a') | strcmp(arg,'ieee-le.l64')
90     ieee='a';
91     elseif strcmp(arg,'s') | strcmp(arg,'ieee-be.l64')
92     ieee='s';
93     elseif strcmp(arg,'rec')
94     userecords=1;
95     else
96     error(['Optional argument ' arg ' is unknown'])
97     end
98     else
99     if userecords==1
100     recnum=arg;
101     elseif isempty(iters)
102     if isnan(arg)
103     iters=scanforfiles(fname);
104     disp([ sprintf('Reading %i time levels:',size(iters,2)) sprintf(' %i',iters) ]);
105     elseif isinf(arg)
106     iters=scanforfiles(fname);
107     if isempty(iters)
108     AA=[];iters=[];return;
109     end
110     disp([ sprintf('Found %i time levels, reading %i',size(iters,2),iters(end)) ]);
111     iters=iters(end);
112     elseif prod(arg>=0) & prod(round(arg)==arg)
113     if arg>=9999999999
114     error(sprintf('Argument %i > 9999999999',arg))
115     end
116     iters=arg;
117     else
118     error(sprintf('Argument %i must be a positive integer',arg))
119     end
120     else
121     error('Multiple iterations should be specified as a vector')
122     end
123     end
124     end
125    
126     if isempty(iters)
127     iters=-1;
128     end
129    
130     % Loop over each iteration
131     for iter=1:size(iters,2);
132     if iters(iter)>=0
133     fname=sprintf('%s.%10.10i',fnamearg,iters(iter));
134     end
135    
136     % Figure out if there is a path in the filename
137     NS=findstr('/',fname);
138     if size(NS)>0
139     Dir=fname(1:NS(end));
140     else
141     Dir='./';
142     end
143    
144     % Match name of all meta-files
145     allfiles=dir( sprintf('%s*.meta',fname) );
146    
147     if size(allfiles,1)==0
148     disp(sprintf('No files match the search: %s*.meta',fname));
149     %allow partial reads% error('No files found.')
150     end
151    
152     % Loop through allfiles
153     for j=1:size(allfiles,1);
154    
155     % Read meta- and data-file
156     [A,N] = localrdmds([Dir allfiles(j).name],ieee,recnum);
157    
158     bdims=N(1,:);
159     r0=N(2,:);
160     rN=N(3,:);
161     ndims=prod(size(bdims));
162     if (ndims == 1)
163     AA(r0(1):rN(1),iter)=A;
164     elseif (ndims == 2)
165     AA(r0(1):rN(1),r0(2):rN(2),iter)=A;
166     elseif (ndims == 3)
167     AA(r0(1):rN(1),r0(2):rN(2),r0(3):rN(3),iter)=A;
168     elseif (ndims == 4)
169     AA(r0(1):rN(1),r0(2):rN(2),r0(3):rN(3),r0(4):rN(4),iter)=A;
170     elseif (ndims == 5)
171     AA(r0(1):rN(1),r0(2):rN(2),r0(3):rN(3),r0(4):rN(4),r0(5):rN(5),iter)=A;
172     else
173     error('Dimension of data set is larger than currently coded. Sorry!')
174     end
175    
176     end % files
177     end % iterations
178    
179     %-------------------------------------------------------------------------------
180    
181     function [A,N] = localrdmds(fname,ieee,recnum)
182    
183     mname=strrep(fname,' ','');
184     dname=strrep(mname,'.meta','.data');
185    
186     % Read and interpret Meta file
187     fid = fopen(mname,'r');
188     if (fid == -1)
189     error(['File' mname ' could not be opened'])
190     end
191    
192     % Scan each line of the Meta file
193     allstr=' ';
194     keepgoing = 1;
195     while keepgoing > 0,
196     line = fgetl(fid);
197     if (line == -1)
198     keepgoing=-1;
199     else
200     % Strip out "(PID.TID *.*)" by finding first ")"
201     %old ind=findstr([line ')'],')'); line=line(ind(1)+1:end);
202     ind=findstr(line,')');
203     if size(ind) ~= 0
204     line=line(ind(1)+1:end);
205     end
206     % Remove comments of form //
207     line=[line ' //']; ind=findstr(line,'//'); line=line(1:ind(1)-1);
208     % Add to total string
209     allstr=[allstr line];
210     end
211     end
212    
213     % Close meta file
214     fclose(fid);
215    
216     % Strip out comments of form /* ... */
217     ind1=findstr(allstr,'/*'); ind2=findstr(allstr,'*/');
218     if size(ind1) ~= size(ind2)
219     error('The /* ... */ comments are not properly paired')
220     end
221     while size(ind1,2) > 0
222     allstr=[allstr(1:ind1(1)-1) allstr(ind2(1)+3:end)];
223     ind1=findstr(allstr,'/*'); ind2=findstr(allstr,'*/');
224     end
225    
226     % This is a kludge to catch whether the meta-file is of the
227     % old or new type. nrecords does not exist in the old type.
228     nrecords = NaN;
229    
230     % Everything in lower case
231     allstr=lower(allstr);
232    
233     % Fix the unfortunate choice of 'format'
234     allstr=strrep(allstr,'format','dataprec');
235    
236     % Evaluate meta information
237     eval(allstr);
238    
239     N=reshape( dimlist , 3 , prod(size(dimlist))/3 );
240     if ~isnan(nrecords) & nrecords > 1 & isempty(recnum)
241     N=[N,[nrecords 1 nrecords]'];
242     elseif ~isempty(recnum) & recnum>nrecords
243     error('Requested record number is higher than the number of available records')
244     end
245    
246     if isempty(recnum)
247     recnum=1;
248     end
249    
250     if isnan(nrecords)
251     % This is the old 'meta' method that used sequential access
252    
253     A=allstr;
254     % Open data file
255     fid=fopen(dname,'r',ieee);
256    
257     % Read record size in bytes
258     recsz=fread(fid,1,'uint32');
259     ldims=N(3,:)-N(2,:)+1;
260     numels=prod(ldims);
261    
262     rat=recsz/numels;
263     if rat == 4
264     A=fread(fid,numels,'real*4');
265     elseif rat == 8
266     A=fread(fid,numels,'real*8');
267     else
268     sprintf(' Implied size in meta-file = %d', numels )
269     sprintf(' Record size in data-file = %d', recsz )
270     error('Ratio between record size and size in meta-file inconsistent')
271     end
272    
273     erecsz=fread(fid,1,'uint32');
274     if erecsz ~= recsz
275     sprintf('WARNING: Record sizes at beginning and end of file are inconsistent')
276     end
277    
278     fclose(fid);
279    
280     A=reshape(A,ldims);
281    
282     else
283     % This is the new MDS format that uses direct access
284    
285     ldims=N(3,:)-N(2,:)+1;
286     for r=1:size(recnum(:),1);
287     if dataprec == 'float32'
288     A(:,r)=myrdda(dname,ldims,recnum(r),'real*4',ieee);
289     elseif dataprec == 'float64'
290     A(:,r)=myrdda(dname,ldims,recnum(r),'real*8',ieee);
291     else
292     error(['Unrecognized format in meta-file = ' format]);
293     end
294     end
295    
296     A=reshape(A,[ldims size(recnum(:),1)]);
297     if size(recnum(:),1)>1
298     N(1,end+1)=size(recnum(:),1);
299     N(2,end)=1;
300     N(3,end)=size(recnum(:),1);
301     end
302    
303     end
304    
305     %-------------------------------------------------------------------------------
306    
307     % result = RDDA( file, dim, irec [options] )
308     %
309     % This routine reads the irec'th record of shape 'dim' from the
310     % direct-access binary file (float or double precision) 'file'.
311     %
312     % Required arguments:
313     %
314     % file - string - name of file to read from
315     % dim - vector - dimensions of the file records and the resulting array
316     % irec - integer - record number in file in which to write data
317     %
318     % Optional arguments (must appear after the required arguments):
319     % prec - string - precision of storage in file. Default = 'real*8'.
320     % ieee - string - IEEE bit-wise representation in file. Default = 'b'.
321     %
322     % 'prec' may take the values:
323     % 'real*4' - floating point, 32 bits.
324     % 'real*8' - floating point, 64 bits - the efault.
325     %
326     % 'ieee' may take values:
327     % 'ieee-be' or 'b' - IEEE floating point with big-endian
328     % byte ordering - the default
329     % 'ieee-le' or 'l' - IEEE floating point with little-endian
330     % byte ordering
331     % 'native' or 'n' - local machine format
332     % 'cray' or 'c' - Cray floating point with big-endian
333     % byte ordering
334     % 'ieee-le.l64' or 'a' - IEEE floating point with little-endian
335     % byte ordering and 64 bit long data type
336     % 'ieee-be.l64' or 's' - IEEE floating point with big-endian byte
337     % ordering and 64 bit long data type.
338     %
339     % eg. T=rdda('t.data',[64 64 32],1);
340     % T=rdda('t.data',[256],4,'real*4');
341     % T=rdda('t.data',[128 64],2,'real*4','b');
342     function [arr] = myrdda(file,N,k,varargin)
343    
344     % Defaults
345     WORDLENGTH=8;
346     rtype='real*8';
347     ieee='b';
348    
349     % Check optional arguments
350     args=char(varargin);
351     while (size(args,1) > 0)
352     if deblank(args(1,:)) == 'real*4'
353     WORDLENGTH=4;
354     rtype='real*4';
355     elseif deblank(args(1,:)) == 'real*8'
356     WORDLENGTH=8;
357     rtype='real*8';
358     elseif deblank(args(1,:)) == 'n' | deblank(args(1,:)) == 'native'
359     ieee='n';
360     elseif deblank(args(1,:)) == 'l' | deblank(args(1,:)) == 'ieee-le'
361     ieee='l';
362     elseif deblank(args(1,:)) == 'b' | deblank(args(1,:)) == 'ieee-be'
363     ieee='b';
364     elseif deblank(args(1,:)) == 'c' | deblank(args(1,:)) == 'cray'
365     ieee='c';
366     elseif deblank(args(1,:)) == 'a' | deblank(args(1,:)) == 'ieee-le.l64'
367     ieee='a';
368     elseif deblank(args(1,:)) == 's' | deblank(args(1,:)) == 'ieee-be.l64'
369     ieee='s';
370     else
371     error(['Optional argument ' args(1,:) ' is unknown'])
372     end
373     args=args(2:end,:);
374     end
375    
376     nnn=prod(N);
377    
378     [fid mess]=fopen(file,'r',ieee);
379     if fid == -1
380     error('Error while opening file:\n%s',mess)
381     end
382     st=fseek(fid,nnn*(k-1)*WORDLENGTH,'bof');
383     if st ~= 0
384     mess=ferror(fid);
385     error('There was an error while positioning the file pointer:\n%s',mess)
386     end
387     [arr count]=fread(fid,nnn,rtype);
388     if count ~= nnn
389     error('Not enough data was available to be read: off EOF?')
390     end
391     st=fclose(fid);
392     %arr=reshape(arr,N);
393    
394     %
395     function [iters] = scanforfiles(fname)
396    
397     iters=[];
398     allfiles=dir([fname '.*.001.001.meta']);
399     if isempty(allfiles)
400     allfiles=dir([fname '.*.meta']);
401     ioff=0;
402     else
403     ioff=8;
404     end
405     for k=1:size(allfiles,1);
406     hh=allfiles(k).name;
407     iters(k)=str2num( hh(end-14-ioff:end-5-ioff) );
408     end
409     iters=sort(iters);

  ViewVC Help
Powered by ViewVC 1.1.22