1 |
function [nIt,rList,tim,vvA,listV,listK]=read_StD(namF,sufx,listV); |
2 |
% [nIt,rList,tim,vvA,listV,listK]=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 |
% listK = list of level numbers |
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.5 2014/05/14 16:53:04 jmc Exp $ |
19 |
% $Name: $ |
20 |
|
21 |
%- Remove insignificant whitespace: |
22 |
%sufx=strtrim(char(sufx)); % <-- only with matlab-7 or more recent |
23 |
sufx=strrep(char(sufx),' ',''); |
24 |
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 |
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 |
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(nIt,1)-tim(1,1))/(nIt-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 |
if nv == 1, listK=nk; else listK=[listK nk]; end |
162 |
end |
163 |
fprintf(' <= end \n'); |
164 |
if length(msgA) ~= 1, fprintf([' not found:',msgA,'\n']) ; end |
165 |
if length(msgB) ~= 1, fprintf([' dim Pb:',msgB,'\n']) ; end |
166 |
if n3d == 0, n3d=1; vvA=zeros(n3d,nIt,nReg,5,nbV); else |
167 |
s3d=max(listK); if n3d > s3d, vvA=vvA(1:s3d,:,:,:,:); end |
168 |
end |
169 |
|
170 |
return |
171 |
|
172 |
function [vvm,nk,msg] = read_StD_1v(namF,nmV,sufx,nReg,nit,dit) |
173 |
msg=' '; nk=0; fprintf([' ',nmV]); |
174 |
namfil=[namF,'_',nmV,'.',sufx]; D=dir(namfil); |
175 |
if size(D,1) == 1, |
176 |
var=load(namfil); |
177 |
nk=1+max(var(:,1)); |
178 |
if size(var,1) == nk*nReg*nit, |
179 |
vv1=reshape(var(:,2:6),[nk*nReg nit 5]); |
180 |
%-- put back Integrated Volume (was cumulated over the time period) in 5: |
181 |
vv1(:,:,5)=vv1(:,:,5).*(ones(nk*nReg,1)*dit); |
182 |
vv1=reshape(vv1,[nk nReg nit 5]); |
183 |
vvm=permute(vv1,[1 3 2 4]); |
184 |
else |
185 |
%fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var)); |
186 |
%fprintf(' ; nk,nReg,nit= %i %i %i\n',nk,nReg,nit); |
187 |
msg=sprintf([' in ',nmV,' : %i <> %ix%ix%i'],size(var,1),nk,nReg,nit); nk=0; |
188 |
vvm=zeros(1,nit,nReg,5); |
189 |
return |
190 |
end |
191 |
else msg=sprintf([msg,namfil,' ']); nk=-1; vvm=zeros(1,nit,nReg,5); |
192 |
end |
193 |
return |