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

Contents of /MITgcm_contrib/jmc_script/read_StD.m

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


Revision 1.3 - (show annotations) (download)
Wed Oct 8 18:18:04 2008 UTC (15 years, 6 months ago) by jmc
Branch: MAIN
Changes since 1.2: +6 -3 lines
update matlab-scripts to plot Stat-Diags output

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

  ViewVC Help
Powered by ViewVC 1.1.22