Warning, /utils/matlab/read_StD.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 9acfcdc5 on 2020-02-12 18:01:28 UTC
9acfcdc58c Jeff*0001 function [nIt,rList,tim,vvA,listV,listK]=read_StD(namF,sufx,listV);
0002 % [nIt,rList,tim,vvA,listV,listK]=read_StD(namF,sufx,listV);
0003 %
0004 % read ASCII stat-Diags output files (after splitted by script extract_StD)
0005 %
0006 % input: listV = list of fileds to read in ; if ='all_flds' => read all fields
0007 % namF = prefix of all file names to read (after extract_StD)
0008 % sufx = suffix of all file names to read (after extract_StD)
0009 % output:
0010 % nIt = number of time reccords
0011 % rList = list of region number
0012 % tim(:,1) = iterations number ; tim(:,2) = time in simulation
0013 % listV = list of fields
0014 % listK = list of level numbers
0015 % vvA = 5 dims output array:
0016 % ( kLev, time_rec, region_rec, [ave,std,min,max,vol], var_rec )
0017
0018 % $Header: /u/gcmpack/MITgcm_contrib/jmc_script/read_StD.m,v 1.6 2019/04/01 22:04:43 jmc Exp $
0019 % $Name: $
0020
0021 %- Remove insignificant whitespace:
0022 %sufx=strtrim(char(sufx)); % <-- only with matlab-7 or more recent
0023 sufx=strrep(char(sufx),' ','');
0024 namfhd=[namF,'_head','.',sufx];
0025 namfil=[namF,'_Iter','.',sufx];
0026 fprintf(['read ',sufx,' :']);
0027 rf=-1; if strcmp(char(listV),'all_flds'), rf=0; end
0028
0029 %- read from Header file:
0030 % frequency ; list of regions ; list of levels ; (+ list of Var, if rf=0)
0031 fid=fopen(namfhd,'r');
0032 D=dir(namfhd);
0033 if size(D,1) == 1,
0034 fprintf(' head');
0035 else error(['file: ',namfhd,' not found']); return; end
0036 flag=1; l=0; k=3;
0037 while flag > 0
0038 tline=fgetl(fid);
0039 if ischar(tline), l=l+1;
0040 %disp(tline);
0041 if tline(2:12) == ' frequency ',
0042 frq=sscanf(tline(17:end),'%f');
0043 k=k-1;
0044 end
0045 if tline(2:10) == ' Regions ',
0046 rList=sscanf(tline(17:end),'%i');
0047 k=k-1;
0048 end
0049 if tline(2:15) == ' Nb of levels ',
0050 kList=sscanf(tline(17:end),'%i');
0051 k=k-1;
0052 end
0053 if rf >=0 & tline(2:10) == ' Fields ',
0054 tmp=[tline(17:end),' ']; i1=0;
0055 for i=1:length(tmp),
0056 if isspace(tmp(i)),
0057 if i1 > 0 & i > i1,
0058 if rf == 0, listV=cellstr(tmp(i1:i-1));
0059 else listV(rf+1)=cellstr(tmp(i1:i-1)); end
0060 i1=i+1; rf=rf+1;
0061 else i1=i+1;
0062 end
0063 end
0064 end
0065 end
0066 if tline(2:15) == ' end of header',
0067 flag=0;
0068 end
0069 else flag=-1; end
0070 end
0071 fclose(fid);
0072 if rf > 0,
0073 %- rename fields (consistent with script "extract_StD"):
0074 for j=1:rf,
0075 var1=char(listV(j));
0076 switch var1
0077 case 'ETAN' , var2='Eta';
0078 case 'ETANSQ' , var2='Et2';
0079 case 'THETA' , var2='T';
0080 case 'SALT' , var2='S';
0081 case 'UVEL' , var2='U';
0082 case 'VVEL' , var2='V';
0083 case 'WVEL' , var2='W';
0084 case 'PHIHYD' , var2='Phi';
0085 case 'UVELSQ' , var2='U2';
0086 case 'VVELSQ' , var2='V2';
0087 case 'WVELSQ' , var2='W2';
0088 case 'THETASQ' , var2='T2';
0089 otherwise var2=var1;
0090 end
0091 listV(j)=cellstr(var2);
0092 if strcmp(var1,var2), fprintf(' %s\n',var2);
0093 else fprintf(' %s --> %s\n',var1,var2); end
0094 end
0095 % listV
0096 end
0097 if flag ~= 0 | k > 0, frq,rList,kList, end
0098 if flag ~= 0, error(['not normal end after reading ',int2str(l),' lines']); end
0099 if k > 0, error(['missing ',int2str(k),' flds (freq,regions,levels)']); end
0100 if flag ~= 0 | k > 0, frq=0; rList=0; kList=0; return; end
0101 nReg=size(rList,1);
0102 nbV=size(listV,2);
0103 %fprintf('\n kList=');fprintf(' %i',kList);fprintf('\n');
0104 %if rf > 0, for l=1:nbV, fprintf([' ',char(listV(l))]); end; fprintf('\n'); end
0105 if ( rf > 0 & nbV < length(kList) ) | nbV > length(kList),
0106 fprintf('\n Warning: nbV= %i but (only) %i in kList\n',nbV,length(kList))
0107 end
0108 %- set Max Nb of levels: n3d
0109 n3d=max(kList); if n3d > 1, n3d=1+n3d ; end
0110
0111 %- load Iter:
0112 D=dir(namfil);
0113 if size(D,1) == 1,
0114 fprintf(' & Iter'); var=load(namfil);
0115 else fprintf(['file: ',namfil,' not found => EXIT \n']); return; end
0116 nIt=size(var,1);
0117
0118 %- initialize output & save Iter:
0119 tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV);
0120 msgA=' '; msgB=' ';
0121 tim(:,1)=var;
0122
0123 dIt=ones(1,nIt);
0124 if nIt > 1, dIt(2:nIt)=tim(2:nIt,1)-tim(1:nIt-1,1); dIt(1)=dIt(2); end
0125 % fprintf(' (dIt: %i %i)',min(dIt(1:nIt)),max(dIt(1:nIt)));
0126 %- take the inverse : will be used to correct the volume
0127 % (frq<0 : no volume correction if snap-shot)
0128 if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end
0129
0130 %- build time:
0131 delT=0;
0132 if nIt > 1, delT=(tim(nIt,1)-tim(1,1))/(nIt-1); delT=max(0,delT); end
0133 if delT > 0,
0134 delT=abs(frq)/delT;
0135 delta=delT-round(delT);
0136 delta=abs(delta*(tim(2,1)-tim(1,1)));
0137 if delta <= 0.5, delT=round(delT); end
0138 end
0139 fprintf(': nReg= %i ; delta.T= %6.1f ; freq= %8.0f',nReg,delT,frq);
0140 if delT > 0, fprintf(' (= %8.3f it)',frq/delT); end
0141 fprintf('\n')
0142
0143 if delT > 0,
0144 tim(:,2)=delT*tim(:,1);
0145 else
0146 tim(:,2)=tim(:,1);
0147 end
0148 fprintf(' + var:'); var=load(namfil);
0149
0150 %---------
0151 for nv=1:nbV,
0152 namV=char(listV(nv));
0153 [vv1,nk,msg]=read_StD_1v(namF,namV,sufx,nReg,nIt,dIt);
0154 if nk > 0, fprintf('_%i',nk); elseif nk < 0, msgA=sprintf([msgA,msg]);
0155 else msgB=sprintf([msgB,msg]) ; end
0156 if nk > n3d,
0157 error([' too many levels(=',int2str(nk),' > ',int2str(n3d), ...
0158 ' reading field: ',namV]);
0159 end
0160 if nk > 0, vvA(1:nk,:,:,:,nv)=vv1; end
0161 if nv == 1, listK=nk; else listK=[listK nk]; end
0162 end
0163 fprintf(' <= end \n');
0164 if length(msgA) ~= 1, fprintf([' not found:',msgA,'\n']) ; end
0165 if length(msgB) ~= 1, fprintf([' dim Pb:',msgB,'\n']) ; end
0166 if n3d == 0, n3d=1; vvA=zeros(n3d,nIt,nReg,5,nbV); else
0167 s3d=max(listK); if n3d > s3d, vvA=vvA(1:s3d,:,:,:,:); end
0168 end
0169
0170 return
0171
0172 function [vvm,nk,msg] = read_StD_1v(namF,nmV,sufx,nReg,nit,dit)
0173 msg=' '; nk=0; fprintf([' ',nmV]);
0174 namfil=[namF,'_',nmV,'.',sufx]; D=dir(namfil);
0175 if size(D,1) == 1,
0176 var=load(namfil);
0177 nk=1+max(var(:,1));
0178 if size(var,1) == nk*nReg*nit,
0179 vv1=reshape(var(:,2:6),[nk*nReg nit 5]);
0180 %-- put back Integrated Volume (was cumulated over the time period) in 5:
0181 vv1(:,:,5)=vv1(:,:,5).*(ones(nk*nReg,1)*dit);
0182 vv1=reshape(vv1,[nk nReg nit 5]);
0183 vvm=permute(vv1,[1 3 2 4]);
0184 else
0185 %fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var));
0186 %fprintf(' ; nk,nReg,nit= %i %i %i\n',nk,nReg,nit);
0187 msg=sprintf([' in ',nmV,' : %i <> %ix%ix%i'],size(var,1),nk,nReg,nit); nk=0;
0188 vvm=zeros(1,nit,nReg,5);
0189 return
0190 end
0191 else msg=sprintf([msg,namfil,' ']); nk=-1; vvm=zeros(1,nit,nReg,5);
0192 end
0193 return
