/[MITgcm]/MITgcm_contrib/jmc_script/read_StD.m
ViewVC logotype

Annotation of /MITgcm_contrib/jmc_script/read_StD.m

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


Revision 1.4 - (hide annotations) (download)
Thu Oct 9 01:09:48 2008 UTC (15 years, 6 months ago) by jmc
Branch: MAIN
Changes since 1.3: +3 -2 lines
compatible with matlab-6

1 jmc 1.3 function [nIt,rList,tim,vvA,listV,kList]=read_StD(namF,sufx,listV);
2     % [nIt,rList,tim,vvA,listV,kList]=read_StD(namF,sufx,listV);
3 jmc 1.1 %
4     % read ASCII stat-Diags output files (after splitted by script extract_StD)
5     %
6     % input: listV = list of fileds to read in ; if ='all_flds' => read all fields
7     % namF = prefix of all file names to read (after extract_StD)
8     % sufx = suffix of all file names to read (after extract_StD)
9     % output:
10     % nIt = number of time reccords
11     % rList = list of region number
12     % tim(:,1) = iterations number ; tim(:,2) = time in simulation
13     % listV = list of fields
14 jmc 1.3 % kList = list of levels number
15 jmc 1.1 % vvA = 5 dims output array:
16     % ( kLev, time_rec, region_rec, [ave,std,min,max,vol], var_rec )
17    
18 jmc 1.4 % $Header: /u/gcmpack/MITgcm_contrib/jmc_script/read_StD.m,v 1.3 2008/10/08 18:18:04 jmc Exp $
19 jmc 1.1 % $Name: $
20    
21 jmc 1.3 %- Remove insignificant whitespace:
22 jmc 1.4 %sufx=strtrim(char(sufx)); % <-- only with matlab-7 or more recent
23     sufx=strrep(char(sufx),' ','');
24 jmc 1.1 namfhd=[namF,'_head','.',sufx];
25     namfil=[namF,'_Iter','.',sufx];
26     fprintf(['read ',sufx,' :']);
27     rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end
28    
29     %- read from Header file:
30     % frequency ; list of regions ; list of levels ; (+ list of Var, if rf=0)
31     fid=fopen(namfhd,'r');
32     D=dir(namfhd);
33     if size(D,1) == 1,
34     fprintf(' head');
35     else error(['file: ',namfhd,' not found']); return; end
36     flag=1; l=0; k=3;
37     while flag > 0
38     tline=fgetl(fid);
39     if ischar(tline), l=l+1;
40     %disp(tline);
41     if tline(2:12) == ' frequency ',
42     frq=sscanf(tline(17:end),'%f');
43     k=k-1;
44     end
45     if tline(2:10) == ' Regions ',
46     rList=sscanf(tline(17:end),'%i');
47     k=k-1;
48     end
49     if tline(2:15) == ' Nb of levels ',
50     kList=sscanf(tline(17:end),'%i');
51     k=k-1;
52     end
53     if rf >=0 & tline(2:10) == ' Fields ',
54     tmp=[tline(17:end),' ']; i1=0;
55     for i=1:length(tmp),
56     if isspace(tmp(i)),
57     if i1 > 0 & i > i1,
58     if rf == 0, listV=cellstr(tmp(i1:i-1));
59     else listV(rf+1)=cellstr(tmp(i1:i-1)); end
60     i1=i+1; rf=rf+1;
61     else i1=i+1;
62     end
63     end
64     end
65     end
66     if tline(2:15) == ' end of header',
67     flag=0;
68     end
69     else flag=-1; end
70     end
71     fclose(fid);
72     if rf > 0,
73     %- rename fields (consistent with script "extract_StD"):
74 jmc 1.2 for j=1:rf,
75     var1=char(listV(j));
76     switch var1
77     case 'ETAN' , var2='Eta';
78     case 'ETANSQ' , var2='Et2';
79     case 'THETA' , var2='T';
80     case 'SALT' , var2='S';
81     case 'UVEL' , var2='U';
82     case 'VVEL' , var2='V';
83     case 'WVEL' , var2='W';
84     case 'PHIHYD' , var2='Phi';
85     case 'UVELSQ' , var2='U2';
86     case 'VVELSQ' , var2='V2';
87     case 'WVELSQ' , var2='W2';
88     case 'THETASQ' , var2='T2';
89     otherwise var2=var1;
90     end
91     listV(j)=cellstr(var2);
92     if strcmp(var1,var2), fprintf(' %s\n',var2);
93     else fprintf(' %s --> %s\n',var1,var2); end
94     end
95     % listV
96 jmc 1.1 end
97     if flag ~= 0 | k > 0, frq,rList,kList, end
98     if flag ~= 0, error(['not normal end after reading ',int2str(l),' lines']); end
99     if k > 0, error(['missing ',int2str(k),' flds (freq,regions,levels)']); end
100     if flag ~= 0 | k > 0, frq=0; rList=0; kList=0; return; end
101     nReg=size(rList,1);
102     nbV=size(listV,2);
103     %fprintf('\n kList=');fprintf(' %i',kList);fprintf('\n');
104     %if rf > 0, for l=1:nbV, fprintf([' ',char(listV(l))]); end; fprintf('\n'); end
105     if ( rf > 0 & nbV < length(kList) ) | nbV > length(kList),
106     fprintf('\n Warning: nbV= %i but (only) %i in kList\n',nbV,length(kList))
107     end
108     %- set Max Nb of levels: n3d
109     n3d=max(kList); if n3d > 1, n3d=1+n3d ; end
110    
111     %- load Iter:
112     D=dir(namfil);
113     if size(D,1) == 1,
114     fprintf(' & Iter'); var=load(namfil);
115     else fprintf(['file: ',namfil,' not found => EXIT \n']); return; end
116     nIt=size(var,1);
117    
118     %- initialize output & save Iter:
119     tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV);
120     msgA=' '; msgB=' ';
121     tim(:,1)=var;
122    
123     dIt=ones(1,nIt);
124     if nIt > 1, dIt(2:nIt)=tim(2:nIt,1)-tim(1:nIt-1,1); dIt(1)=dIt(2); end
125     % fprintf(' (dIt: %i %i)',min(dIt(1:nIt)),max(dIt(1:nIt)));
126     %- take the inverse : will be used to correct the volume
127     % (frq<0 : no volume correction if snap-shot)
128     if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end
129    
130     %- build time:
131     delT=0;
132     if nIt > 1, delT=tim(2,1)-tim(1,1); delT=max(0,delT); end
133     if delT > 0,
134     delT=abs(frq)/delT;
135     delta=delT-round(delT);
136     delta=abs(delta*(tim(2,1)-tim(1,1)));
137     if delta <= 0.5, delT=round(delT); end
138     end
139     fprintf(': nReg= %i ; delta.T= %6.1f ; freq= %8.0f',nReg,delT,frq);
140     if delT > 0, fprintf(' (= %8.3f it)',frq/delT); end
141     fprintf('\n')
142    
143     if delT > 0,
144     tim(:,2)=delT*tim(:,1);
145     else
146     tim(:,2)=tim(:,1);
147     end
148     fprintf(' + var:'); var=load(namfil);
149    
150     %---------
151     for nv=1:nbV,
152     namV=char(listV(nv));
153     [vv1,nk,msg]=read_StD_1v(namF,namV,sufx,nReg,nIt,dIt);
154     if nk > 0, fprintf('_%i',nk); elseif nk < 0, msgA=sprintf([msgA,msg]);
155     else msgB=sprintf([msgB,msg]) ; end
156     if nk > n3d,
157     error([' too many levels(=',int2str(nk),' > ',int2str(n3d), ...
158     ' reading field: ',namV]);
159     end
160     if nk > 0, vvA(1:nk,:,:,:,nv)=vv1; end
161     end
162     fprintf(' <= end \n');
163     if length(msgA) ~= 1, fprintf([' not found:',msgA,'\n']) ; end
164     if length(msgB) ~= 1, fprintf([' dim Pb:',msgB,'\n']) ; end
165     if n3d == 0, n3d=1; vvA=zeros(n3d,nIt,nReg,5,nbV); end
166    
167     return
168    
169     function [vvm,nk,msg] = read_StD_1v(namF,nmV,sufx,nReg,nit,dit)
170     msg=' '; nk=0; fprintf([' ',nmV]);
171     namfil=[namF,'_',nmV,'.',sufx]; D=dir(namfil);
172     if size(D,1) == 1,
173     var=load(namfil);
174     nk=1+max(var(:,1));
175     if size(var,1) == nk*nReg*nit,
176     vv1=reshape(var(:,2:6),[nk*nReg nit 5]);
177     %-- put back Integrated Volume (was cumulated over the time period) in 5:
178     vv1(:,:,5)=vv1(:,:,5).*(ones(nk*nReg,1)*dit);
179     vv1=reshape(vv1,[nk nReg nit 5]);
180     vvm=permute(vv1,[1 3 2 4]);
181     else
182     %fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var));
183     %fprintf(' ; nk,nReg,nit= %i %i %i\n',nk,nReg,nit);
184     msg=sprintf([' in ',nmV,' : %i <> %ix%ix%i'],size(var,1),nk,nReg,nit); nk=0;
185     vvm=zeros(1,nit,nReg,5);
186     return
187     end
188     else msg=sprintf([msg,namfil,' ']); nk=-1; vvm=zeros(1,nit,nReg,5);
189     end
190     return

  ViewVC Help
Powered by ViewVC 1.1.22