| 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 |
jmc |
1.2 |
% $Header: /u/gcmpack/MITgcm_contrib/jmc_script/read_StD.m,v 1.1 2007/03/01 04:35:35 jmc Exp $ |
| 18 |
jmc |
1.1 |
% $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 |
jmc |
1.2 |
for j=1:rf, |
| 71 |
|
|
var1=char(listV(j)); |
| 72 |
|
|
switch var1 |
| 73 |
|
|
case 'ETAN' , var2='Eta'; |
| 74 |
|
|
case 'ETANSQ' , var2='Et2'; |
| 75 |
|
|
case 'THETA' , var2='T'; |
| 76 |
|
|
case 'SALT' , var2='S'; |
| 77 |
|
|
case 'UVEL' , var2='U'; |
| 78 |
|
|
case 'VVEL' , var2='V'; |
| 79 |
|
|
case 'WVEL' , var2='W'; |
| 80 |
|
|
case 'PHIHYD' , var2='Phi'; |
| 81 |
|
|
case 'UVELSQ' , var2='U2'; |
| 82 |
|
|
case 'VVELSQ' , var2='V2'; |
| 83 |
|
|
case 'WVELSQ' , var2='W2'; |
| 84 |
|
|
case 'THETASQ' , var2='T2'; |
| 85 |
|
|
otherwise var2=var1; |
| 86 |
|
|
end |
| 87 |
|
|
listV(j)=cellstr(var2); |
| 88 |
|
|
if strcmp(var1,var2), fprintf(' %s\n',var2); |
| 89 |
|
|
else fprintf(' %s --> %s\n',var1,var2); end |
| 90 |
|
|
end |
| 91 |
|
|
% listV |
| 92 |
jmc |
1.1 |
end |
| 93 |
|
|
if flag ~= 0 | k > 0, frq,rList,kList, end |
| 94 |
|
|
if flag ~= 0, error(['not normal end after reading ',int2str(l),' lines']); end |
| 95 |
|
|
if k > 0, error(['missing ',int2str(k),' flds (freq,regions,levels)']); end |
| 96 |
|
|
if flag ~= 0 | k > 0, frq=0; rList=0; kList=0; return; end |
| 97 |
|
|
nReg=size(rList,1); |
| 98 |
|
|
nbV=size(listV,2); |
| 99 |
|
|
%fprintf('\n kList=');fprintf(' %i',kList);fprintf('\n'); |
| 100 |
|
|
%if rf > 0, for l=1:nbV, fprintf([' ',char(listV(l))]); end; fprintf('\n'); end |
| 101 |
|
|
if ( rf > 0 & nbV < length(kList) ) | nbV > length(kList), |
| 102 |
|
|
fprintf('\n Warning: nbV= %i but (only) %i in kList\n',nbV,length(kList)) |
| 103 |
|
|
end |
| 104 |
|
|
%- set Max Nb of levels: n3d |
| 105 |
|
|
n3d=max(kList); if n3d > 1, n3d=1+n3d ; end |
| 106 |
|
|
|
| 107 |
|
|
%- load Iter: |
| 108 |
|
|
D=dir(namfil); |
| 109 |
|
|
if size(D,1) == 1, |
| 110 |
|
|
fprintf(' & Iter'); var=load(namfil); |
| 111 |
|
|
else fprintf(['file: ',namfil,' not found => EXIT \n']); return; end |
| 112 |
|
|
nIt=size(var,1); |
| 113 |
|
|
|
| 114 |
|
|
%- initialize output & save Iter: |
| 115 |
|
|
tim=zeros(nIt,2); vvA=zeros(n3d,nIt,nReg,5,nbV); |
| 116 |
|
|
msgA=' '; msgB=' '; |
| 117 |
|
|
tim(:,1)=var; |
| 118 |
|
|
|
| 119 |
|
|
dIt=ones(1,nIt); |
| 120 |
|
|
if nIt > 1, dIt(2:nIt)=tim(2:nIt,1)-tim(1:nIt-1,1); dIt(1)=dIt(2); end |
| 121 |
|
|
% fprintf(' (dIt: %i %i)',min(dIt(1:nIt)),max(dIt(1:nIt))); |
| 122 |
|
|
%- take the inverse : will be used to correct the volume |
| 123 |
|
|
% (frq<0 : no volume correction if snap-shot) |
| 124 |
|
|
if frq > 0, dIt=1./max(1,dIt); else dIt=ones(1,nIt); end |
| 125 |
|
|
|
| 126 |
|
|
%- build time: |
| 127 |
|
|
delT=0; |
| 128 |
|
|
if nIt > 1, delT=tim(2,1)-tim(1,1); delT=max(0,delT); end |
| 129 |
|
|
if delT > 0, |
| 130 |
|
|
delT=abs(frq)/delT; |
| 131 |
|
|
delta=delT-round(delT); |
| 132 |
|
|
delta=abs(delta*(tim(2,1)-tim(1,1))); |
| 133 |
|
|
if delta <= 0.5, delT=round(delT); end |
| 134 |
|
|
end |
| 135 |
|
|
fprintf(': nReg= %i ; delta.T= %6.1f ; freq= %8.0f',nReg,delT,frq); |
| 136 |
|
|
if delT > 0, fprintf(' (= %8.3f it)',frq/delT); end |
| 137 |
|
|
fprintf('\n') |
| 138 |
|
|
|
| 139 |
|
|
if delT > 0, |
| 140 |
|
|
tim(:,2)=delT*tim(:,1); |
| 141 |
|
|
else |
| 142 |
|
|
tim(:,2)=tim(:,1); |
| 143 |
|
|
end |
| 144 |
|
|
fprintf(' + var:'); var=load(namfil); |
| 145 |
|
|
|
| 146 |
|
|
%--------- |
| 147 |
|
|
for nv=1:nbV, |
| 148 |
|
|
namV=char(listV(nv)); |
| 149 |
|
|
[vv1,nk,msg]=read_StD_1v(namF,namV,sufx,nReg,nIt,dIt); |
| 150 |
|
|
if nk > 0, fprintf('_%i',nk); elseif nk < 0, msgA=sprintf([msgA,msg]); |
| 151 |
|
|
else msgB=sprintf([msgB,msg]) ; end |
| 152 |
|
|
if nk > n3d, |
| 153 |
|
|
error([' too many levels(=',int2str(nk),' > ',int2str(n3d), ... |
| 154 |
|
|
' reading field: ',namV]); |
| 155 |
|
|
end |
| 156 |
|
|
if nk > 0, vvA(1:nk,:,:,:,nv)=vv1; end |
| 157 |
|
|
end |
| 158 |
|
|
fprintf(' <= end \n'); |
| 159 |
|
|
if length(msgA) ~= 1, fprintf([' not found:',msgA,'\n']) ; end |
| 160 |
|
|
if length(msgB) ~= 1, fprintf([' dim Pb:',msgB,'\n']) ; end |
| 161 |
|
|
if n3d == 0, n3d=1; vvA=zeros(n3d,nIt,nReg,5,nbV); end |
| 162 |
|
|
|
| 163 |
|
|
return |
| 164 |
|
|
|
| 165 |
|
|
function [vvm,nk,msg] = read_StD_1v(namF,nmV,sufx,nReg,nit,dit) |
| 166 |
|
|
msg=' '; nk=0; fprintf([' ',nmV]); |
| 167 |
|
|
namfil=[namF,'_',nmV,'.',sufx]; D=dir(namfil); |
| 168 |
|
|
if size(D,1) == 1, |
| 169 |
|
|
var=load(namfil); |
| 170 |
|
|
nk=1+max(var(:,1)); |
| 171 |
|
|
if size(var,1) == nk*nReg*nit, |
| 172 |
|
|
vv1=reshape(var(:,2:6),[nk*nReg nit 5]); |
| 173 |
|
|
%-- put back Integrated Volume (was cumulated over the time period) in 5: |
| 174 |
|
|
vv1(:,:,5)=vv1(:,:,5).*(ones(nk*nReg,1)*dit); |
| 175 |
|
|
vv1=reshape(vv1,[nk nReg nit 5]); |
| 176 |
|
|
vvm=permute(vv1,[1 3 2 4]); |
| 177 |
|
|
else |
| 178 |
|
|
%fprintf(['\n ERROR in Dim, var=',nmV,' :']); fprintf(' %i',size(var)); |
| 179 |
|
|
%fprintf(' ; nk,nReg,nit= %i %i %i\n',nk,nReg,nit); |
| 180 |
|
|
msg=sprintf([' in ',nmV,' : %i <> %ix%ix%i'],size(var,1),nk,nReg,nit); nk=0; |
| 181 |
|
|
vvm=zeros(1,nit,nReg,5); |
| 182 |
|
|
return |
| 183 |
|
|
end |
| 184 |
|
|
else msg=sprintf([msg,namfil,' ']); nk=-1; vvm=zeros(1,nit,nReg,5); |
| 185 |
|
|
end |
| 186 |
|
|
return |
Parent Directory
| 
Revision Log
| 
Revision Graph