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)

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

if lChunk==0; return; end;

dirModel=[dirModel '/'];

if isempty(dirMat); dirMat=[dirModel 'mat/']; end;
if isempty(dir(dirMat)); eval(['!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/' ]};

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

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