--- MITgcm_contrib/jmc_script/read_StD.m	2007/03/01 04:35:35	1.1
+++ MITgcm_contrib/jmc_script/read_StD.m	2019/04/01 22:04:43	1.6
@@ -1,40 +1,44 @@
-function [nIt,rList,tim,vvA,listV]=read_StD(namF,sufx,listV);
-% [nIt,rList,tim,vvA,listV]=read_StD(namF,sufx,listV);
+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: 
+% 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: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/jmc_script/read_StD.m,v 1.1 2007/03/01 04:35:35 jmc Exp $
+% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/jmc_script/read_StD.m,v 1.6 2019/04/01 22:04:43 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: 
+%- 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; 
+flag=1; l=0; k=3;
 while flag > 0
  tline=fgetl(fid);
- if ischar(tline), l=l+1; 
+ if ischar(tline), l=l+1;
    %disp(tline);
-    if tline(2:12) == ' frequency ', 
+    if tline(2:12) == ' frequency ',
      frq=sscanf(tline(17:end),'%f');
      k=k-1;
     end
@@ -42,24 +46,24 @@
      rList=sscanf(tline(17:end),'%i');
      k=k-1;
     end
-    if tline(2:15) == ' Nb of levels ', 
+    if tline(2:15) == ' Nb of levels ',
      kList=sscanf(tline(17:end),'%i');
      k=k-1;
     end
-    if rf >=0 & tline(2:10) == ' Fields  ', 
+    if rf >=0 & tline(2:10) == ' Fields  ',
       tmp=[tline(17:end),' ']; i1=0;
       for i=1:length(tmp),
-        if isspace(tmp(i)), 
+        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; 
+          else i1=i+1;
           end
         end
       end
     end
-    if tline(2:15) == ' end of header', 
+    if tline(2:15) == ' end of header',
      flag=0;
     end
  else flag=-1; end
@@ -67,11 +71,28 @@
 fclose(fid);
 if rf > 0,
 %- rename fields (consistent with script "extract_StD"):
-   listV=strrep(listV,'ETAN', 'Eta'); listV=strrep(listV,'ETANSQ','Et2');
-   listV=strrep(listV,'THETA' ,'T');  listV=strrep(listV,'SALT'  ,'S');
-   listV=strrep(listV,'UVEL'  ,'U');  listV=strrep(listV,'VVEL'  ,'V');
-   listV=strrep(listV,'WVEL'  ,'W');  listV=strrep(listV,'PHIHYD','Phi');
-   listV=strrep(listV,'UVELSQ','U2'); listV=strrep(listV,'VVELSQ','V2');
+   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
@@ -97,8 +118,8 @@
 %- initialize output & save Iter:
 tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV);
 msgA=' '; msgB=' ';
-tim(:,1)=var; 
- 
+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)));
@@ -106,9 +127,9 @@
 %            (frq<0 : no volume correction if snap-shot)
 if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end
 
-%- build time: 
+%- build time:
 delT=0;
-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
 if delT > 0,
   delT=abs(frq)/delT;
   delta=delT-round(delT);
@@ -137,11 +158,14 @@
            ' 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); 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
 
@@ -158,7 +182,7 @@
      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(['\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);