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

Diff of /MITgcm_contrib/jmc_script/read_StD.m

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

revision 1.1 by jmc, Thu Mar 1 04:35:35 2007 UTC revision 1.6 by jmc, Mon Apr 1 22:04:43 2019 UTC
# Line 1  Line 1 
1  function [nIt,rList,tim,vvA,listV]=read_StD(namF,sufx,listV);  function [nIt,rList,tim,vvA,listV,listK]=read_StD(namF,sufx,listV);
2  % [nIt,rList,tim,vvA,listV]=read_StD(namF,sufx,listV);  % [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)  % 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  % 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)  %        namF  = prefix of all file names to read (after extract_StD)
8  %        sufx  = suffix of all file names to read (after extract_StD)  %        sufx  = suffix of all file names to read (after extract_StD)
9  % output:  % output:
10  %  nIt      = number of time reccords  %  nIt      = number of time reccords
11  %  rList    = list of region number  %  rList    = list of region number
12  %  tim(:,1) = iterations number ; tim(:,2) = time in simulation  %  tim(:,1) = iterations number ; tim(:,2) = time in simulation
13  %  listV    = list of fields  %  listV    = list of fields
14    %  listK    = list of level numbers
15  %  vvA      = 5 dims output array:  %  vvA      = 5 dims output array:
16  %           ( kLev, time_rec, region_rec, [ave,std,min,max,vol], var_rec )  %           ( kLev, time_rec, region_rec, [ave,std,min,max,vol], var_rec )
17    
18  % $Header$  % $Header$
19  % $Name$  % $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];  namfhd=[namF,'_head','.',sufx];
25  namfil=[namF,'_Iter','.',sufx];  namfil=[namF,'_Iter','.',sufx];
26   fprintf(['read ',sufx,' :']);   fprintf(['read ',sufx,' :']);
27  rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end  rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end
28    
29  %- read from Header file:  %- read from Header file:
30  %   frequency ; list of regions ; list of levels ; (+ list of Var, if rf=0)  %   frequency ; list of regions ; list of levels ; (+ list of Var, if rf=0)
31  fid=fopen(namfhd,'r');  fid=fopen(namfhd,'r');
32  D=dir(namfhd);  D=dir(namfhd);
33  if size(D,1) == 1,  if size(D,1) == 1,
34   fprintf(' head');   fprintf(' head');
35  else error(['file: ',namfhd,' not found']); return; end  else error(['file: ',namfhd,' not found']); return; end
36  flag=1; l=0; k=3;  flag=1; l=0; k=3;
37  while flag > 0  while flag > 0
38   tline=fgetl(fid);   tline=fgetl(fid);
39   if ischar(tline), l=l+1;   if ischar(tline), l=l+1;
40     %disp(tline);     %disp(tline);
41      if tline(2:12) == ' frequency ',      if tline(2:12) == ' frequency ',
42       frq=sscanf(tline(17:end),'%f');       frq=sscanf(tline(17:end),'%f');
43       k=k-1;       k=k-1;
44      end      end
# Line 42  while flag > 0 Line 46  while flag > 0
46       rList=sscanf(tline(17:end),'%i');       rList=sscanf(tline(17:end),'%i');
47       k=k-1;       k=k-1;
48      end      end
49      if tline(2:15) == ' Nb of levels ',      if tline(2:15) == ' Nb of levels ',
50       kList=sscanf(tline(17:end),'%i');       kList=sscanf(tline(17:end),'%i');
51       k=k-1;       k=k-1;
52      end      end
53      if rf >=0 & tline(2:10) == ' Fields  ',      if rf >=0 & tline(2:10) == ' Fields  ',
54        tmp=[tline(17:end),' ']; i1=0;        tmp=[tline(17:end),' ']; i1=0;
55        for i=1:length(tmp),        for i=1:length(tmp),
56          if isspace(tmp(i)),          if isspace(tmp(i)),
57            if i1 > 0 & i > i1,            if i1 > 0 & i > i1,
58              if rf == 0, listV=cellstr(tmp(i1:i-1));              if rf == 0, listV=cellstr(tmp(i1:i-1));
59              else listV(rf+1)=cellstr(tmp(i1:i-1)); end              else listV(rf+1)=cellstr(tmp(i1:i-1)); end
60              i1=i+1; rf=rf+1;              i1=i+1; rf=rf+1;
61            else i1=i+1;            else i1=i+1;
62            end            end
63          end          end
64        end        end
65      end      end
66      if tline(2:15) == ' end of header',      if tline(2:15) == ' end of header',
67       flag=0;       flag=0;
68      end      end
69   else flag=-1; end   else flag=-1; end
# Line 67  end Line 71  end
71  fclose(fid);  fclose(fid);
72  if rf > 0,  if rf > 0,
73  %- rename fields (consistent with script "extract_StD"):  %- rename fields (consistent with script "extract_StD"):
74     listV=strrep(listV,'ETAN', 'Eta'); listV=strrep(listV,'ETANSQ','Et2');     for j=1:rf,
75     listV=strrep(listV,'THETA' ,'T');  listV=strrep(listV,'SALT'  ,'S');       var1=char(listV(j));
76     listV=strrep(listV,'UVEL'  ,'U');  listV=strrep(listV,'VVEL'  ,'V');       switch var1
77     listV=strrep(listV,'WVEL'  ,'W');  listV=strrep(listV,'PHIHYD','Phi');       case 'ETAN'    ,  var2='Eta';
78     listV=strrep(listV,'UVELSQ','U2'); listV=strrep(listV,'VVELSQ','V2');       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  end
97  if flag ~= 0 | k > 0, frq,rList,kList, end  if flag ~= 0 | k > 0, frq,rList,kList, end
98  if flag ~= 0, error(['not normal end after reading ',int2str(l),' lines']); end  if flag ~= 0, error(['not normal end after reading ',int2str(l),' lines']); end
# Line 97  nIt=size(var,1); Line 118  nIt=size(var,1);
118  %- initialize output & save Iter:  %- initialize output & save Iter:
119  tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV);  tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV);
120  msgA=' '; msgB=' ';  msgA=' '; msgB=' ';
121  tim(:,1)=var;  tim(:,1)=var;
122    
123  dIt=ones(1,nIt);  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  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)));  % fprintf(' (dIt: %i %i)',min(dIt(1:nIt)),max(dIt(1:nIt)));
# Line 106  if nIt > 1, dIt(2:nIt)=tim(2:nIt,1)-tim( Line 127  if nIt > 1, dIt(2:nIt)=tim(2:nIt,1)-tim(
127  %            (frq<0 : no volume correction if snap-shot)  %            (frq<0 : no volume correction if snap-shot)
128  if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end  if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end
129    
130  %- build time:  %- build time:
131  delT=0;  delT=0;
132  if nIt > 1, delT=tim(2,1)-tim(1,1); delT=max(0,delT); end  if nIt > 1, delT=(tim(nIt,1)-tim(1,1))/(nIt-1); delT=max(0,delT); end
133  if delT > 0,  if delT > 0,
134    delT=abs(frq)/delT;    delT=abs(frq)/delT;
135    delta=delT-round(delT);    delta=delT-round(delT);
# Line 137  for nv=1:nbV, Line 158  for nv=1:nbV,
158             ' reading field: ',namV]);             ' reading field: ',namV]);
159    end    end
160    if nk > 0, vvA(1:nk,:,:,:,nv)=vv1; end    if nk > 0, vvA(1:nk,:,:,:,nv)=vv1; end
161      if nv == 1, listK=nk; else listK=[listK nk]; end
162  end  end
163  fprintf(' <= end \n');  fprintf(' <= end \n');
164  if length(msgA) ~= 1, fprintf(['  not found:',msgA,'\n']) ; end  if length(msgA) ~= 1, fprintf(['  not found:',msgA,'\n']) ; end
165  if length(msgB) ~= 1, fprintf(['  dim Pb:',msgB,'\n']) ; end  if length(msgB) ~= 1, fprintf(['  dim Pb:',msgB,'\n']) ; end
166  if n3d == 0, n3d=1; vvA=zeros(n3d,nIt,nReg,5,nbV); end  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  return
171    
# Line 158  function [vvm,nk,msg] = read_StD_1v(namF Line 182  function [vvm,nk,msg] = read_StD_1v(namF
182       vv1=reshape(vv1,[nk nReg nit 5]);       vv1=reshape(vv1,[nk nReg nit 5]);
183       vvm=permute(vv1,[1 3 2 4]);       vvm=permute(vv1,[1 3 2 4]);
184     else     else
185      %fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var));      %fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var));
186      %fprintf(' ; nk,nReg,nit= %i %i %i\n',nk,nReg,nit);      %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;       msg=sprintf([' in ',nmV,' : %i <> %ix%ix%i'],size(var,1),nk,nReg,nit); nk=0;
188       vvm=zeros(1,nit,nReg,5);       vvm=zeros(1,nit,nReg,5);
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22