/[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.2 - (show annotations) (download)
Thu May 29 23:47:13 2008 UTC (15 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.1: +23 -6 lines
improve scripts to read & plot Stat-Diags (ascii) output files

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

  ViewVC Help
Powered by ViewVC 1.1.22