matlab_class/gcmfaces_diags/diags_select.m

function []=basic_diags_ecco(dirModel,dirMat,setDiags,lChunk,iChunk,varargin);
%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
%              iChunk states which part of the records should be
%                     computed (i0+1:i0+lChunk where i0~lTot/lChunk)

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

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;

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


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