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;

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