matlab_class/gcmfaces_diags/diags_select.m

function []=diags_select(dirModel,dirMat,setDiags,lChunk,listChunk);
%object:       compute a set of basic diagnostics
%input:        dirModel is the model run directory
%              dirMat is the directory where diagnozed .mat files will be saved
%                     -> set it to '' to use the default [dirModel 'mat/']
%              setDiags is the choice of diagnostics
%                       'A') trasnports (see diags_set_A.m)
%                       'B') air-sea fluxes
%                       'C') state variables
%                       'D') global and regional budgets
%                       'LAYERS') T/S/RHO coordinate overturns
%                       'user') is a placeholder for user diags additions
%              lChunk states how many records are computed at once
%              listChunk states which parts of the records should be
%                     computed (i0+1:i0+lChunk where i0~lTot/lChunk)

gcmfaces_global;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%create dirModel/mat if needed:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if lChunk==0; return; end;

dirModel=[dirModel '/'];

if isempty(dirMat); dirMat=[dirModel 'mat/']; end;
if isempty(dir(dirMat)); mkdir(dirMat); end;

%list of subdirectories where to scan for output
listSubdirs={[dirMat 'BUDG/' ],[dirModel 'diags/BUDG/' ],[dirModel 'diags/OTHER/' ],...
    [dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/'],[dirModel 'diags/' ]};

if isfield(myenv,'nctiles');
  use_nctiles=myenv.nctiles;
else;
  use_nctiles=0;
end;

setDiagsParams=[];
if iscell(setDiags);
    setDiagsParams={setDiags{2:end}};
    setDiags=setDiags{1};
end;

%%%%%%%%%%%%%%%%%%%%%%
%load grid and params:
%%%%%%%%%%%%%%%%%%%%%%

gcmfaces_global; global myparms;
test0=~isempty(dir([dirMat 'diags_grid_parms.mat']));
test1=~isempty(dir([dirMat 'basic_diags_ecco_mygrid.mat']));
if test0;%load mygrid and myparms from file
    eval(['load ' dirMat 'diags_grid_parms.mat;']);
elseif test1;%for backward compatibility
    eval(['load ' dirMat 'basic_diags_ecco_mygrid.mat;']);
else;
    error('could not find diags_grid_parms.mat');
end;

%in case mygrid.memoryLimit=1, load the stuff that was not saved to diags_grid_parms.mat
if mygrid.memoryLimit==1;
        list0={'hFacS','hFacW'};
        for iFld=1:length(list0);
          eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([mygrid.dirGrid ''' list0{iFld} '*'']);']);
        end;
        %
        mygrid.hFacCsurf=mygrid.hFacC;
        for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
        %
        mskC=mygrid.hFacC; mskC(mskC==0)=NaN; mskC(mskC>0)=1; mygrid.mskC=mskC;
        mskW=mygrid.hFacW; mskW(mskW==0)=NaN; mskW(mskW>0)=1; mygrid.mskW=mskW;
        mskS=mygrid.hFacS; mskS(mskS==0)=NaN; mskS(mskS>0)=1; mygrid.mskS=mskS;
        %
        gcmfaces_lines_zonal;
        mygrid.LATS=[mygrid.LATS_MASKS.lat]';
        [lonPairs,latPairs,names]=line_greatC_TUV_MASKS_v4;
        gcmfaces_lines_transp(lonPairs,latPairs,names);
end;

%%%%%%%%%%%%%%
%define diags:
%%%%%%%%%%%%%%

userStep=1;
if ~isempty(which(['diags_set_' setDiags]));
    eval(['diags_set_' setDiags]);
else;
    diags_set_user;
end;

%reformat listDiags to cell array:
jj=strfind(listDiags,' '); jj=[0 jj length(listDiags)+1];
for ii=1:length(jj)-1;
    tmp1=listDiags(jj(ii)+1:jj(ii+1)-1);
    if ii==1; listDiags2={tmp1}; else; listDiags2{ii}=tmp1; end;
end;
listDiags=listDiags2; clear listDiags2;

%%%%%%%%%%%%%%
%detect files:
%%%%%%%%%%%%%%

userStep=2;
if ~isempty(which(['diags_set_' setDiags]));
    eval(['diags_set_' setDiags]);
else;
    diags_set_user;
end;

%remove blanks in pkg/diagnostics flds name
listFlds=deblank(listFlds); 
listFldsNames=deblank(listFldsNames);

%set the list of diags times
[listTimes]=diags_list_times(listSubdirs,listFiles);

%then scan directories (listSubdirs) for files (listFiles)
if ~use_nctiles;
for ii=1:length(listFiles);
    tmp0='';
    for jj=1:length(listSubdirs);
        tmp1=listFiles{ii}; tmp2=listSubdirs{jj};
        if ~isempty(dir([tmp2 '/' tmp1 '*meta']));
            tmp0=[tmp2 '/' tmp1]; listFiles{ii}=tmp0;
        end;
    end;
    if isempty(tmp0); fprintf([' not found: ' tmp1 '\n']); listFiles{ii}=''; end;
end;
end;

if use_nctiles;
for ii=1:length(listFlds);
tmp0=sum(strcmp(myenv.nctileslist,listFlds{ii}));
if tmp0==0; error([' not found: ' tmp1 '\n']); end;
if tmp0<1; error([' found several : ' tmp1 '\n']); end;
end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%loop over chunks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isdir([dirMat 'diags_set_' setDiags]);
  dirMat=[dirMat 'diags_set_' setDiags '/'];
end;

for iChunk=listChunk;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%load records for computation:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lTot=length(listTimes);
i0=min(lChunk*(iChunk-1),lTot);
i1=min(i0+lChunk,lTot);

tic;

if ~use_nctiles;
tt=listTimes([i0+1:i1])';
for iFile=1:length(listFiles);
    fileFld=listFiles{iFile};
    if ~isempty(fileFld);
        rdmds2workspace_list(fileFld,tt,listFlds);
        for ii=1:length(meta.fldList);
            tmp1=deblank(meta.fldList{ii});
            jj=find(strcmp(tmp1,listFlds));
            if ~isempty(jj); eval(['all_' listFldsNames{jj} '=' tmp1 ';']); end;
            eval(['clear ' tmp1]);
        end;
    end;
end;
end;

if use_nctiles;
tt=[i0+1:i1]';
di=myenv.nctilesdir;
for iFld=1:length(listFlds);
  nm=listFlds{iFld};
  fileIn=sprintf('%s/%s/%s',di,nm,nm);
  fld=read_nctiles(fileIn,nm,tt);
  eval(['all_' nm '=fld;']);
end;
end;

fprintf([num2str(i1-i0) ' records loaded in ' num2str(toc) '\n']);

%check that all fields were found, and otherwise reduce listFlds
listFldsIsPresent=zeros(1,length(listFldsNames));
for iFld=1:length(listFldsNames);
    fldName=['all_' listFldsNames{iFld}];
    if isempty(who(fldName));
        fprintf([fldName ' is missing \n']);
    else;
        listFldsIsPresent(iFld)=1;
        %    fprintf([fldName ' was found \n']);
    end;
end;
listFldsNames={listFldsNames{find(listFldsIsPresent)}};
listFlds={listFlds{find(listFldsIsPresent)}};

%%%%%%%%%%%%%%%%%%%%
%do the computation:
%%%%%%%%%%%%%%%%%%%%
for ii=i0+1:i1;
    
    tic;
    tt=listTimes(ii);
    
    %get data from buffer
    for jj=1:length(listFldsNames);
        if lChunk==1|length(listTimes)==1;
            eval([listFldsNames{jj} '=all_' listFldsNames{jj} ';']);
        else;
            eval(['tmp1=size(all_' listFldsNames{jj} '{1});']);
            if length(tmp1)==4; eval([listFldsNames{jj} '=all_' listFldsNames{jj} '(:,:,:,ii-i0);']);
            else; eval([listFldsNames{jj} '=all_' listFldsNames{jj} '(:,:,ii-i0);']);
            end;
        end;
    end;
    fprintf([num2str(ii-i0) '/' num2str(i1-i0) ' started \n']);

    %generic output file name (that will potentially be overriden)
    tmp1=setDiags;
    fileMat=['diags_set_' tmp1 '_' num2str(tt) '.mat'];
        
    %the actual computations
    userStep=3;
    if ~isempty(which(['diags_set_' setDiags]));
        eval(['diags_set_' setDiags]);
    else;
        diags_set_user;
    end;
   
    fprintf([num2str(ii) '/' num2str(length(listTimes)) ' done in ' num2str(toc) '\n']);
 
    %package results
    onediag=[];
    onediag.listTimes=myparms.yearFirst(1)+tt*myparms.timeStep/86400/365.25;
    onediag.listSteps=tt;
    for jj=1:length(listDiags);
        eval(['onediag.' listDiags{jj} '=' listDiags{jj} ';']);
    end;

    %write to disk 
    diags_store_to_mat(dirMat,fileMat,onediag);
    
end;%for ii=i0+1:i1;

end;%for iChunk=listChunk;


1	gforget	1.3	function []=diags_select(dirModel,dirMat,setDiags,lChunk,listChunk);
2	gforget	1.1	%object: compute a set of basic diagnostics
3			%input: dirModel is the model run directory
4			% dirMat is the directory where diagnozed .mat files will be saved
5			% -> set it to '' to use the default [dirModel 'mat/']
6			% setDiags is the choice of diagnostics
7			% 'A') trasnports (see diags_set_A.m)
8			% 'B') air-sea fluxes
9			% 'C') state variables
10			% 'D') global and regional budgets
11			% 'LAYERS') T/S/RHO coordinate overturns
12			% 'user') is a placeholder for user diags additions
13			% lChunk states how many records are computed at once
14	gforget	1.3	% listChunk states which parts of the records should be
15	gforget	1.1	% computed (i0+1:i0+lChunk where i0~lTot/lChunk)
16
17	gforget	1.5	gcmfaces_global;
18
19	gforget	1.1	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20			%create dirModel/mat if needed:
21			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22
23	gforget	1.3	if lChunk==0; return; end;
24
25	gforget	1.1	dirModel=[dirModel '/'];
26
27			if isempty(dirMat); dirMat=[dirModel 'mat/']; end;
28	gforget	1.6	if isempty(dir(dirMat)); mkdir(dirMat); end;
29	gforget	1.1
30			%list of subdirectories where to scan for output
31			listSubdirs={[dirMat 'BUDG/' ],[dirModel 'diags/BUDG/' ],[dirModel 'diags/OTHER/' ],...
32			[dirModel 'diags/STATE/' ],[dirModel 'diags/TRSP/'],[dirModel 'diags/' ]};
33
34	gforget	1.5	if isfield(myenv,'nctiles');
35			use_nctiles=myenv.nctiles;
36			else;
37			use_nctiles=0;
38			end;
39
40	gforget	1.1	setDiagsParams=[];
41			if iscell(setDiags);
42			setDiagsParams={setDiags{2:end}};
43			setDiags=setDiags{1};
44			end;
45
46			%%%%%%%%%%%%%%%%%%%%%%
47			%load grid and params:
48			%%%%%%%%%%%%%%%%%%%%%%
49
50			gcmfaces_global; global myparms;
51			test0=~isempty(dir([dirMat 'diags_grid_parms.mat']));
52			test1=~isempty(dir([dirMat 'basic_diags_ecco_mygrid.mat']));
53			if test0;%load mygrid and myparms from file
54			eval(['load ' dirMat 'diags_grid_parms.mat;']);
55			elseif test1;%for backward compatibility
56			eval(['load ' dirMat 'basic_diags_ecco_mygrid.mat;']);
57			else;
58			error('could not find diags_grid_parms.mat');
59			end;
60
61	gforget	1.2	%in case mygrid.memoryLimit=1, load the stuff that was not saved to diags_grid_parms.mat
62			if mygrid.memoryLimit==1;
63			list0={'hFacS','hFacW'};
64			for iFld=1:length(list0);
65			eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([mygrid.dirGrid ''' list0{iFld} '*'']);']);
66			end;
67			%
68			mygrid.hFacCsurf=mygrid.hFacC;
69			for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
70			%
71			mskC=mygrid.hFacC; mskC(mskC==0)=NaN; mskC(mskC>0)=1; mygrid.mskC=mskC;
72			mskW=mygrid.hFacW; mskW(mskW==0)=NaN; mskW(mskW>0)=1; mygrid.mskW=mskW;
73			mskS=mygrid.hFacS; mskS(mskS==0)=NaN; mskS(mskS>0)=1; mygrid.mskS=mskS;
74			%
75			gcmfaces_lines_zonal;
76			mygrid.LATS=[mygrid.LATS_MASKS.lat]';
77			[lonPairs,latPairs,names]=line_greatC_TUV_MASKS_v4;
78			gcmfaces_lines_transp(lonPairs,latPairs,names);
79			end;
80
81	gforget	1.1	%%%%%%%%%%%%%%
82			%define diags:
83			%%%%%%%%%%%%%%
84
85			userStep=1;
86			if ~isempty(which(['diags_set_' setDiags]));
87			eval(['diags_set_' setDiags]);
88			else;
89			diags_set_user;
90			end;
91
92			%reformat listDiags to cell array:
93			jj=strfind(listDiags,' '); jj=[0 jj length(listDiags)+1];
94			for ii=1:length(jj)-1;
95			tmp1=listDiags(jj(ii)+1:jj(ii+1)-1);
96			if ii==1; listDiags2={tmp1}; else; listDiags2{ii}=tmp1; end;
97			end;
98			listDiags=listDiags2; clear listDiags2;
99
100			%%%%%%%%%%%%%%
101			%detect files:
102			%%%%%%%%%%%%%%
103
104			userStep=2;
105			if ~isempty(which(['diags_set_' setDiags]));
106			eval(['diags_set_' setDiags]);
107			else;
108			diags_set_user;
109			end;
110
111			%remove blanks in pkg/diagnostics flds name
112			listFlds=deblank(listFlds);
113			listFldsNames=deblank(listFldsNames);
114
115			%set the list of diags times
116			[listTimes]=diags_list_times(listSubdirs,listFiles);
117
118			%then scan directories (listSubdirs) for files (listFiles)
119	gforget	1.5	if ~use_nctiles;
120	gforget	1.1	for ii=1:length(listFiles);
121			tmp0='';
122			for jj=1:length(listSubdirs);
123			tmp1=listFiles{ii}; tmp2=listSubdirs{jj};
124			if ~isempty(dir([tmp2 '/' tmp1 '*meta']));
125			tmp0=[tmp2 '/' tmp1]; listFiles{ii}=tmp0;
126			end;
127			end;
128			if isempty(tmp0); fprintf([' not found: ' tmp1 '\n']); listFiles{ii}=''; end;
129			end;
130	gforget	1.5	end;
131
132			if use_nctiles;
133			for ii=1:length(listFlds);
134			tmp0=sum(strcmp(myenv.nctileslist,listFlds{ii}));
135			if tmp0==0; error([' not found: ' tmp1 '\n']); end;
136			if tmp0<1; error([' found several : ' tmp1 '\n']); end;
137			end;
138			end;
139	gforget	1.1
140			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141	gforget	1.3	%loop over chunks
142			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
143
144	gforget	1.4	if isdir([dirMat 'diags_set_' setDiags]);
145			dirMat=[dirMat 'diags_set_' setDiags '/'];
146			end;
147
148	gforget	1.3	for iChunk=listChunk;
149
150			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151	gforget	1.1	%load records for computation:
152			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153
154			lTot=length(listTimes);
155			i0=min(lChunk*(iChunk-1),lTot);
156			i1=min(i0+lChunk,lTot);
157
158			tic;
159	gforget	1.5
160			if ~use_nctiles;
161	gforget	1.1	tt=listTimes([i0+1:i1])';
162			for iFile=1:length(listFiles);
163			fileFld=listFiles{iFile};
164			if ~isempty(fileFld);
165			rdmds2workspace_list(fileFld,tt,listFlds);
166			for ii=1:length(meta.fldList);
167			tmp1=deblank(meta.fldList{ii});
168			jj=find(strcmp(tmp1,listFlds));
169			if ~isempty(jj); eval(['all_' listFldsNames{jj} '=' tmp1 ';']); end;
170			eval(['clear ' tmp1]);
171			end;
172			end;
173			end;
174	gforget	1.5	end;
175
176			if use_nctiles;
177			tt=[i0+1:i1]';
178			di=myenv.nctilesdir;
179			for iFld=1:length(listFlds);
180			nm=listFlds{iFld};
181			fileIn=sprintf('%s/%s/%s',di,nm,nm);
182			fld=read_nctiles(fileIn,nm,tt);
183			eval(['all_' nm '=fld;']);
184			end;
185			end;
186
187	gforget	1.1	fprintf([num2str(i1-i0) ' records loaded in ' num2str(toc) '\n']);
188
189			%check that all fields were found, and otherwise reduce listFlds
190			listFldsIsPresent=zeros(1,length(listFldsNames));
191			for iFld=1:length(listFldsNames);
192			fldName=['all_' listFldsNames{iFld}];
193			if isempty(who(fldName));
194			fprintf([fldName ' is missing \n']);
195			else;
196			listFldsIsPresent(iFld)=1;
197			% fprintf([fldName ' was found \n']);
198			end;
199			end;
200			listFldsNames={listFldsNames{find(listFldsIsPresent)}};
201			listFlds={listFlds{find(listFldsIsPresent)}};
202
203			%%%%%%%%%%%%%%%%%%%%
204			%do the computation:
205			%%%%%%%%%%%%%%%%%%%%
206			for ii=i0+1:i1;
207
208			tic;
209			tt=listTimes(ii);
210
211			%get data from buffer
212			for jj=1:length(listFldsNames);
213			if lChunk==1\|length(listTimes)==1;
214			eval([listFldsNames{jj} '=all_' listFldsNames{jj} ';']);
215			else;
216			eval(['tmp1=size(all_' listFldsNames{jj} '{1});']);
217			if length(tmp1)==4; eval([listFldsNames{jj} '=all_' listFldsNames{jj} '(:,:,:,ii-i0);']);
218			else; eval([listFldsNames{jj} '=all_' listFldsNames{jj} '(:,:,ii-i0);']);
219			end;
220			end;
221			end;
222			fprintf([num2str(ii-i0) '/' num2str(i1-i0) ' started \n']);
223
224			%generic output file name (that will potentially be overriden)
225			tmp1=setDiags;
226			fileMat=['diags_set_' tmp1 '_' num2str(tt) '.mat'];
227
228			%the actual computations
229			userStep=3;
230			if ~isempty(which(['diags_set_' setDiags]));
231			eval(['diags_set_' setDiags]);
232			else;
233			diags_set_user;
234			end;
235	gforget	1.2
236			fprintf([num2str(ii) '/' num2str(length(listTimes)) ' done in ' num2str(toc) '\n']);
237
238	gforget	1.1	%package results
239			onediag=[];
240			onediag.listTimes=myparms.yearFirst(1)+tt*myparms.timeStep/86400/365.25;
241			onediag.listSteps=tt;
242			for jj=1:length(listDiags);
243			eval(['onediag.' listDiags{jj} '=' listDiags{jj} ';']);
244			end;
245
246			%write to disk
247			diags_store_to_mat(dirMat,fileMat,onediag);
248
249	gforget	1.3	end;%for ii=i0+1:i1;
250
251			end;%for iChunk=listChunk;
252	gforget	1.1
253