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

Annotation of /MITgcm_contrib/jmc_script/read_StD.m

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


Revision 1.1 - (hide annotations) (download)
Thu Mar 1 04:35:35 2007 UTC (17 years, 2 months ago) by jmc
Branch: MAIN
add a matlab script to load ASCII atat-diags output files

1 jmc 1.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