matlab_class/gcmfaces_diags/diags_select.m

function [fldsIn]=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)

%to do: add second optional output argument for fldsOut

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;

test0=~isempty(dir([dirMat 'diags_select_budget_list.mat']));
if test0; 
    eval(['load ' dirMat 'diags_select_budget_list.mat;']);
else;
    budget_list=1;
end;
myparms.budgetList=budget_list;

%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;

listSubdirs={[dirModel 'diags/'],listSubdirs{:}};

%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 isempty(strfind(dirMat,['diags_set_' setDiags]));
 dirMat=[dirMat 'diags_set_' setDiags filesep];
 if isempty(dir(dirMat)); mkdir(dirMat); end;
end;

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