MITgcm_contrib/jmc_script/read_StD.m

function [nIt,rList,tim,vvA,listV,listK]=read_StD(namF,sufx,listV);
% [nIt,rList,tim,vvA,listV,listK]=read_StD(namF,sufx,listV);
%
% read ASCII stat-Diags output files (after splitted by script extract_StD)
%
% input: listV = list of fileds to read in ; if ='all_flds' => read all fields
%        namF  = prefix of all file names to read (after extract_StD)
%        sufx  = suffix of all file names to read (after extract_StD)
% output:
%  nIt      = number of time reccords
%  rList    = list of region number
%  tim(:,1) = iterations number ; tim(:,2) = time in simulation
%  listV    = list of fields
%  listK    = list of level numbers
%  vvA      = 5 dims output array:
%           ( kLev, time_rec, region_rec, [ave,std,min,max,vol], var_rec )

% $Header: /u/gcmpack/MITgcm_contrib/jmc_script/read_StD.m,v 1.5 2014/05/14 16:53:04 jmc Exp $
% $Name:  $

%- Remove insignificant whitespace:
%sufx=strtrim(char(sufx)); % <-- only with matlab-7 or more recent
sufx=strrep(char(sufx),' ','');
namfhd=[namF,'_head','.',sufx];
namfil=[namF,'_Iter','.',sufx];
 fprintf(['read ',sufx,' :']);
rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end

%- read from Header file:
%   frequency ; list of regions ; list of levels ; (+ list of Var, if rf=0)
fid=fopen(namfhd,'r');
D=dir(namfhd);
if size(D,1) == 1,
 fprintf(' head');
else error(['file: ',namfhd,' not found']); return; end
flag=1; l=0; k=3;
while flag > 0
 tline=fgetl(fid);
 if ischar(tline), l=l+1;
   %disp(tline);
    if tline(2:12) == ' frequency ',
     frq=sscanf(tline(17:end),'%f');
     k=k-1;
    end
    if tline(2:10) == ' Regions ',
     rList=sscanf(tline(17:end),'%i');
     k=k-1;
    end
    if tline(2:15) == ' Nb of levels ',
     kList=sscanf(tline(17:end),'%i');
     k=k-1;
    end
    if rf >=0 & tline(2:10) == ' Fields  ',
      tmp=[tline(17:end),' ']; i1=0;
      for i=1:length(tmp),
        if isspace(tmp(i)),
          if i1 > 0 & i > i1,
            if rf == 0, listV=cellstr(tmp(i1:i-1));
            else listV(rf+1)=cellstr(tmp(i1:i-1)); end
            i1=i+1; rf=rf+1;
          else i1=i+1;
          end
        end
      end
    end
    if tline(2:15) == ' end of header',
     flag=0;
    end
 else flag=-1; end
end
fclose(fid);
if rf > 0,
%- rename fields (consistent with script "extract_StD"):
   for j=1:rf,
     var1=char(listV(j));
     switch var1
     case 'ETAN'    ,  var2='Eta';
     case 'ETANSQ'  ,  var2='Et2';
     case 'THETA'   ,  var2='T';
     case 'SALT'    ,  var2='S';
     case 'UVEL'    ,  var2='U';
     case 'VVEL'    ,  var2='V';
     case 'WVEL'    ,  var2='W';
     case 'PHIHYD'  ,  var2='Phi';
     case 'UVELSQ'  ,  var2='U2';
     case 'VVELSQ'  ,  var2='V2';
     case 'WVELSQ'  ,  var2='W2';
     case 'THETASQ' ,  var2='T2';
     otherwise     var2=var1;
     end
     listV(j)=cellstr(var2);
     if strcmp(var1,var2), fprintf(' %s\n',var2);
          else fprintf(' %s --> %s\n',var1,var2); end
   end
%  listV
end
if flag ~= 0 | k > 0, frq,rList,kList, end
if flag ~= 0, error(['not normal end after reading ',int2str(l),' lines']); end
if k > 0, error(['missing ',int2str(k),' flds (freq,regions,levels)']); end
if flag ~= 0 | k > 0, frq=0; rList=0; kList=0; return; end
nReg=size(rList,1);
nbV=size(listV,2);
%fprintf('\n kList=');fprintf(' %i',kList);fprintf('\n');
%if rf > 0, for l=1:nbV, fprintf([' ',char(listV(l))]); end; fprintf('\n'); end
if ( rf > 0 & nbV < length(kList) ) | nbV > length(kList),
 fprintf('\n Warning: nbV= %i but (only) %i in kList\n',nbV,length(kList))
end
%- set Max Nb of levels: n3d
n3d=max(kList); if n3d > 1, n3d=1+n3d ; end

%- load Iter:
D=dir(namfil);
if size(D,1) == 1,
 fprintf(' & Iter'); var=load(namfil);
else fprintf(['file: ',namfil,' not found => EXIT \n']); return; end
nIt=size(var,1);

%- initialize output & save Iter:
tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV);
msgA=' '; msgB=' ';
tim(:,1)=var;

dIt=ones(1,nIt);
if nIt > 1, dIt(2:nIt)=tim(2:nIt,1)-tim(1:nIt-1,1); dIt(1)=dIt(2); end
% fprintf(' (dIt: %i %i)',min(dIt(1:nIt)),max(dIt(1:nIt)));
%- take the inverse : will be used to correct the volume
%            (frq<0 : no volume correction if snap-shot)
if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end

%- build time:
delT=0;
if nIt > 1, delT=(tim(nIt,1)-tim(1,1))/(nIt-1); delT=max(0,delT); end
if delT > 0,
  delT=abs(frq)/delT;
  delta=delT-round(delT);
  delta=abs(delta*(tim(2,1)-tim(1,1)));
  if delta <= 0.5, delT=round(delT); end
end
fprintf(': nReg= %i ; delta.T= %6.1f ; freq= %8.0f',nReg,delT,frq);
if delT > 0, fprintf(' (= %8.3f it)',frq/delT); end
fprintf('\n')

if delT > 0,
 tim(:,2)=delT*tim(:,1);
else
 tim(:,2)=tim(:,1);
end
 fprintf(' + var:'); var=load(namfil);

%---------
for nv=1:nbV,
  namV=char(listV(nv));
 [vv1,nk,msg]=read_StD_1v(namF,namV,sufx,nReg,nIt,dIt);
  if nk > 0, fprintf('_%i',nk); elseif nk < 0, msgA=sprintf([msgA,msg]);
  else msgB=sprintf([msgB,msg]) ; end
  if nk > n3d,
    error([' too many levels(=',int2str(nk),' > ',int2str(n3d), ...
           ' reading field: ',namV]);
  end
  if nk > 0, vvA(1:nk,:,:,:,nv)=vv1; end
  if nv == 1, listK=nk; else listK=[listK nk]; end
end
fprintf(' <= end \n');
if length(msgA) ~= 1, fprintf(['  not found:',msgA,'\n']) ; end
if length(msgB) ~= 1, fprintf(['  dim Pb:',msgB,'\n']) ; end
if n3d == 0, n3d=1; vvA=zeros(n3d,nIt,nReg,5,nbV); else
  s3d=max(listK); if n3d > s3d, vvA=vvA(1:s3d,:,:,:,:); end
end

return

function [vvm,nk,msg] = read_StD_1v(namF,nmV,sufx,nReg,nit,dit)
 msg=' '; nk=0; fprintf([' ',nmV]);
 namfil=[namF,'_',nmV,'.',sufx]; D=dir(namfil);
 if size(D,1) == 1,
   var=load(namfil);
   nk=1+max(var(:,1));
   if size(var,1) == nk*nReg*nit,
     vv1=reshape(var(:,2:6),[nk*nReg nit 5]);
%--  put back Integrated Volume (was cumulated over the time period) in 5:
     vv1(:,:,5)=vv1(:,:,5).*(ones(nk*nReg,1)*dit);
     vv1=reshape(vv1,[nk nReg nit 5]);
     vvm=permute(vv1,[1 3 2 4]);
   else
    %fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var));
    %fprintf(' ; nk,nReg,nit= %i %i %i\n',nk,nReg,nit);
     msg=sprintf([' in ',nmV,' : %i <> %ix%ix%i'],size(var,1),nk,nReg,nit); nk=0;
     vvm=zeros(1,nit,nReg,5);
     return
   end
 else msg=sprintf([msg,namfil,' ']); nk=-1; vvm=zeros(1,nit,nReg,5);
 end
return
1	jmc	1.6	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	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	jmc	1.6	% output:
10	jmc	1.1	% 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.6	% listK = list of level numbers
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.6	% $Header: /u/gcmpack/MITgcm_contrib/jmc_script/read_StD.m,v 1.5 2014/05/14 16:53: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	jmc	1.6	%- read from Header file:
30	jmc	1.1	% 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	jmc	1.6	flag=1; l=0; k=3;
37	jmc	1.1	while flag > 0
38			tline=fgetl(fid);
39	jmc	1.6	if ischar(tline), l=l+1;
40	jmc	1.1	%disp(tline);
41	jmc	1.6	if tline(2:12) == ' frequency ',
42	jmc	1.1	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	jmc	1.6	if tline(2:15) == ' Nb of levels ',
50	jmc	1.1	kList=sscanf(tline(17:end),'%i');
51			k=k-1;
52			end
53	jmc	1.6	if rf >=0 & tline(2:10) == ' Fields ',
54	jmc	1.1	tmp=[tline(17:end),' ']; i1=0;
55			for i=1:length(tmp),
56	jmc	1.6	if isspace(tmp(i)),
57	jmc	1.1	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	jmc	1.6	else i1=i+1;
62	jmc	1.1	end
63			end
64			end
65			end
66	jmc	1.6	if tline(2:15) == ' end of header',
67	jmc	1.1	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	jmc	1.6	listV(j)=cellstr(var2);
92	jmc	1.2	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	jmc	1.6	tim(:,1)=var;
122
123	jmc	1.1	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	jmc	1.6	%- build time:
131	jmc	1.1	delT=0;
132	jmc	1.5	if nIt > 1, delT=(tim(nIt,1)-tim(1,1))/(nIt-1); delT=max(0,delT); end
133	jmc	1.1	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	jmc	1.6	if nv == 1, listK=nk; else listK=[listK nk]; end
162	jmc	1.1	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	jmc	1.6	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	jmc	1.1
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) == nknRegnit,
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(nknReg,1)*dit);
182			vv1=reshape(vv1,[nk nReg nit 5]);
183			vvm=permute(vv1,[1 3 2 4]);
184			else
185	jmc	1.6	%fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var));
186	jmc	1.1	%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