/[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.1 - (show annotations) (download)
Thu Mar 1 04:35:35 2007 UTC (17 years, 1 month ago) by jmc
Branch: MAIN
add a matlab script to load ASCII atat-diags 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: $
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 listV=strrep(listV,'ETAN', 'Eta'); listV=strrep(listV,'ETANSQ','Et2');
71 listV=strrep(listV,'THETA' ,'T'); listV=strrep(listV,'SALT' ,'S');
72 listV=strrep(listV,'UVEL' ,'U'); listV=strrep(listV,'VVEL' ,'V');
73 listV=strrep(listV,'WVEL' ,'W'); listV=strrep(listV,'PHIHYD','Phi');
74 listV=strrep(listV,'UVELSQ','U2'); listV=strrep(listV,'VVELSQ','V2');
75 end
76 if flag ~= 0 | k > 0, frq,rList,kList, end
77 if flag ~= 0, error(['not normal end after reading ',int2str(l),' lines']); end
78 if k > 0, error(['missing ',int2str(k),' flds (freq,regions,levels)']); end
79 if flag ~= 0 | k > 0, frq=0; rList=0; kList=0; return; end
80 nReg=size(rList,1);
81 nbV=size(listV,2);
82 %fprintf('\n kList=');fprintf(' %i',kList);fprintf('\n');
83 %if rf > 0, for l=1:nbV, fprintf([' ',char(listV(l))]); end; fprintf('\n'); end
84 if ( rf > 0 & nbV < length(kList) ) | nbV > length(kList),
85 fprintf('\n Warning: nbV= %i but (only) %i in kList\n',nbV,length(kList))
86 end
87 %- set Max Nb of levels: n3d
88 n3d=max(kList); if n3d > 1, n3d=1+n3d ; end
89
90 %- load Iter:
91 D=dir(namfil);
92 if size(D,1) == 1,
93 fprintf(' & Iter'); var=load(namfil);
94 else fprintf(['file: ',namfil,' not found => EXIT \n']); return; end
95 nIt=size(var,1);
96
97 %- initialize output & save Iter:
98 tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV);
99 msgA=' '; msgB=' ';
100 tim(:,1)=var;
101
102 dIt=ones(1,nIt);
103 if nIt > 1, dIt(2:nIt)=tim(2:nIt,1)-tim(1:nIt-1,1); dIt(1)=dIt(2); end
104 % fprintf(' (dIt: %i %i)',min(dIt(1:nIt)),max(dIt(1:nIt)));
105 %- take the inverse : will be used to correct the volume
106 % (frq<0 : no volume correction if snap-shot)
107 if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end
108
109 %- build time:
110 delT=0;
111 if nIt > 1, delT=tim(2,1)-tim(1,1); delT=max(0,delT); end
112 if delT > 0,
113 delT=abs(frq)/delT;
114 delta=delT-round(delT);
115 delta=abs(delta*(tim(2,1)-tim(1,1)));
116 if delta <= 0.5, delT=round(delT); end
117 end
118 fprintf(': nReg= %i ; delta.T= %6.1f ; freq= %8.0f',nReg,delT,frq);
119 if delT > 0, fprintf(' (= %8.3f it)',frq/delT); end
120 fprintf('\n')
121
122 if delT > 0,
123 tim(:,2)=delT*tim(:,1);
124 else
125 tim(:,2)=tim(:,1);
126 end
127 fprintf(' + var:'); var=load(namfil);
128
129 %---------
130 for nv=1:nbV,
131 namV=char(listV(nv));
132 [vv1,nk,msg]=read_StD_1v(namF,namV,sufx,nReg,nIt,dIt);
133 if nk > 0, fprintf('_%i',nk); elseif nk < 0, msgA=sprintf([msgA,msg]);
134 else msgB=sprintf([msgB,msg]) ; end
135 if nk > n3d,
136 error([' too many levels(=',int2str(nk),' > ',int2str(n3d), ...
137 ' reading field: ',namV]);
138 end
139 if nk > 0, vvA(1:nk,:,:,:,nv)=vv1; end
140 end
141 fprintf(' <= end \n');
142 if length(msgA) ~= 1, fprintf([' not found:',msgA,'\n']) ; end
143 if length(msgB) ~= 1, fprintf([' dim Pb:',msgB,'\n']) ; end
144 if n3d == 0, n3d=1; vvA=zeros(n3d,nIt,nReg,5,nbV); end
145
146 return
147
148 function [vvm,nk,msg] = read_StD_1v(namF,nmV,sufx,nReg,nit,dit)
149 msg=' '; nk=0; fprintf([' ',nmV]);
150 namfil=[namF,'_',nmV,'.',sufx]; D=dir(namfil);
151 if size(D,1) == 1,
152 var=load(namfil);
153 nk=1+max(var(:,1));
154 if size(var,1) == nk*nReg*nit,
155 vv1=reshape(var(:,2:6),[nk*nReg nit 5]);
156 %-- put back Integrated Volume (was cumulated over the time period) in 5:
157 vv1(:,:,5)=vv1(:,:,5).*(ones(nk*nReg,1)*dit);
158 vv1=reshape(vv1,[nk nReg nit 5]);
159 vvm=permute(vv1,[1 3 2 4]);
160 else
161 %fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var));
162 %fprintf(' ; nk,nReg,nit= %i %i %i\n',nk,nReg,nit);
163 msg=sprintf([' in ',nmV,' : %i <> %ix%ix%i'],size(var,1),nk,nReg,nit); nk=0;
164 vvm=zeros(1,nit,nReg,5);
165 return
166 end
167 else msg=sprintf([msg,namfil,' ']); nk=-1; vvm=zeros(1,nit,nReg,5);
168 end
169 return

  ViewVC Help
Powered by ViewVC 1.1.22