matlab_class/gcmfaces_diags/diags_display.m

function []=diags_display(dirMat,setDiags,dirTex,nameTex);
%object:        display a set of diagnostics
%input:         dirMat is the directory where diagnozed .mat files are.
%                  dirMat is usually specified as a chararcter string, but
%                  user can also specify dirMat as {dirMat,dirMatRef} 
%                  in order to plot the dirMat-dirMatRef anlomalies.
%               setDiags is the choice of diagnostics set
%                       'A') trasnports
%                       'B') air-sea fluxes
%                       'C') state variables
%                       'D') global and regional budgets
%optional:      dirTex is the directory where tex and figures files are created
%                 (if not specified then display all results to screen instead)
%               nameTex is the tex file name (default : 'myPlots')
%note: specifying dirTex / nameTex requires that the tex file has been 
%      initialized (see e.g. diags_driver_tex.m). Each figure will then 
%      be save as an eps file, added to the tex file, and closed.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%determine input/output params:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%directory names:
if iscell(dirMat); dirMatRef=dirMat{2}; dirMat=dirMat{1}; end;
dirMat=[dirMat '/'];
if isempty(who('dirMatRef')); dirMatRef=''; 
elseif ~isempty(dirMatRef); dirMatRef=[dirMatRef '/']; 
end;
if isempty(who('dirTex')); dirTex=''; else; dirTex=[dirTex '/']; end;
if isempty(who('nameTex')); nameTex='myPlots'; end;

%determine if and where to create tex and figures files
if ~ischar(dirTex); error('mis-specified dirTex'); end;
if isempty(dirTex);
    addToTex=0;
else;
    addToTex=1; fileTex=[dirTex nameTex '.tex'];
end;

%determined where to display anomalies between runs
doAnomalies=~isempty(dirMatRef);

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

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

gcmfaces_global; global myparms;
test1=~isempty(dir([dirMat 'basic_diags_ecco_mygrid.mat']));
test2=~isempty(dir([dirMat 'diags_grid_parms.mat']));
if ~test1&~test2;
  error('missing diags_grid_parms.mat')
elseif test2;
  nameGrid='diags_grid_parms.mat';
  suffDiag='diags_set_';
  budgetList='diags_select_budget_list.mat';
else;
  nameGrid='basic_diags_ecco_mygrid.mat';
  suffDiag='basic_diags_ecco_';
  budgetList='basic_diags_ecco_budget_list.mat';
end;

%here we always reload the grid from dirMat to make sure the same one is used throughout
eval(['load ' dirMat nameGrid ';']);

%backward compatibility:
if ~isfield(mygrid,'memoryLimit'); mygrid.memoryLimit=0; end;
if ~isfield(mygrid,'ioSize'); mygrid.ioSize=0; 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:
%%%%%%%%%%%%%%

if strcmp(nameGrid,'diags_grid_parms.mat');
  %load listDiags and get fileMat (if applies)
  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;
end;

%set fileMat to default name if needed
if isempty(who('fileMat')); fileMat=[suffDiag setDiags]; end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%are there results to display?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

test1=isempty(dir([dirMat '/' fileMat '_*.mat']));
tmp1=[dirMat '/' fileMat '_*.mat']; tmp2=strfind(tmp1,'_'); tmp1(tmp2)=' ';
if test1&addToTex;
    eval(['load ' dirTex '/myPlots.mat;']);
    write2tex(fileTex,3,{'[ ',mySection,' ]'},1);
    write2tex(fileTex,3,{'abort : did not find any'},1);
    write2tex(fileTex,3,{tmp1},1);
    write2tex(fileTex,3,{'results files to display'},1);
    return;
elseif test1;
    fprintf(['\n  abort : did not find any \n    ' tmp1 '\n    results files to display\n\n']);
    return;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%
%load pre-computed diags:
%%%%%%%%%%%%%%%%%%%%%%%%%

tic;

diagsWereLoaded=0;

%specific load sequence (if any, diagsWereLoaded will be set to 1)
if strcmp(nameGrid,'diags_grid_parms.mat');
  userStep=0;
  if ~isempty(which(['diags_set_' setDiags]));
    eval(['diags_set_' setDiags]);
  else;
    diags_set_user;
  end;
end;

%generic load (if no specific one, or )
if ~diagsWereLoaded;
  alldiag=alldiag_load(dirMat,[fileMat '_*.mat']);
end;

toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%load ref and compute anomalies:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if doAnomalies;
    %     scaleAnom=10;
    scaleAnom=1;

    %store:
    dirMatBak=dirMat; dirMat=dirMatRef; alldiagBak=alldiag;

    %load reference solution:
    diagsWereLoaded=0;

    %specific load sequence (if any, diagsWereLoaded will be set to 1)
    if strcmp(nameGrid,'diags_grid_parms.mat');
      userStep=0;
      if ~isempty(which(['diags_set_' setDiags]));
        eval(['diags_set_' setDiags]);
      else;
        diags_set_user;
      end;
    end;

    %generic load (if no specific one, or )
    if ~diagsWereLoaded;
      alldiag=alldiag_load(dirMat,[fileMat '_*.mat']);
    end;

    %restore:
    dirMat=dirMatBak; alldiagRef=alldiag; alldiag=alldiagBak;

    %compute anomalies:
    for ii=1:length(alldiag.listDiags);
        tmp0=alldiag.listDiags{ii};
        if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps')
        if isfield(alldiagRef,tmp0);%compute difference
            tmp1=getfield(alldiag,tmp0)-getfield(alldiagRef,tmp0);
            alldiag=setfield(alldiag,tmp0,tmp1);
        else;%cannot compute diff -> set to NaN
            tmp1=NaN*getfield(alldiag,tmp0);
            alldiag=setfield(alldiag,tmp0,tmp1);
        end;
        end;
    end;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%determine time parameters for plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%number of months for runmean
if myparms.diagsAreMonthly;
  myNmean=12; myNmeanTxt=[' -- ' num2str(myNmean*2) ' months low pass filtered'];
else;
  myNmean=0; myNmeanTxt='';
end;

%first and last year of time average
myYmean=myparms.yearInAve;
if myYmean(1)==myYmean(2); myYmeanTxt=[num2str(myYmean(1)) ' '];
else; myYmeanTxt=[num2str(myYmean(1)) '-' num2str(myYmean(2)) ' '];
end;

%to compute time mean/std on the fly
tt=[myparms.recInAve(1):min(myparms.recInAve(2),length(alldiag.listTimes))];
TT=alldiag.listTimes(tt);
nt=length(TT);

%if only one record, swicth off time series and variance plots
multiTimes=1*(myparms.recInAve(2)>myparms.recInAve(1));

%if llc90 we can plot overturn etc per basin
multiBasins=(sum([90 1170]~=mygrid.ioSize)==0);

%%%%%%%%%%%%%%%
%display diags:
%%%%%%%%%%%%%%%

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


1	gforget	1.1	function []=diags_display(dirMat,setDiags,dirTex,nameTex);
2			%object: display a set of diagnostics
3			%input: dirMat is the directory where diagnozed .mat files are.
4			% dirMat is usually specified as a chararcter string, but
5			% user can also specify dirMat as {dirMat,dirMatRef}
6			% in order to plot the dirMat-dirMatRef anlomalies.
7			% setDiags is the choice of diagnostics set
8			% 'A') trasnports
9			% 'B') air-sea fluxes
10			% 'C') state variables
11			% 'D') global and regional budgets
12			%optional: dirTex is the directory where tex and figures files are created
13			% (if not specified then display all results to screen instead)
14			% nameTex is the tex file name (default : 'myPlots')
15			%note: specifying dirTex / nameTex requires that the tex file has been
16			% initialized (see e.g. diags_driver_tex.m). Each figure will then
17			% be save as an eps file, added to the tex file, and closed.
18
19			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20			%determine input/output params:
21			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22
23			%directory names:
24			if iscell(dirMat); dirMatRef=dirMat{2}; dirMat=dirMat{1}; end;
25			dirMat=[dirMat '/'];
26			if isempty(who('dirMatRef')); dirMatRef='';
27			elseif ~isempty(dirMatRef); dirMatRef=[dirMatRef '/'];
28			end;
29			if isempty(who('dirTex')); dirTex=''; else; dirTex=[dirTex '/']; end;
30			if isempty(who('nameTex')); nameTex='myPlots'; end;
31
32			%determine if and where to create tex and figures files
33			if ~ischar(dirTex); error('mis-specified dirTex'); end;
34			if isempty(dirTex);
35			addToTex=0;
36			else;
37			addToTex=1; fileTex=[dirTex nameTex '.tex'];
38			end;
39
40			%determined where to display anomalies between runs
41			doAnomalies=~isempty(dirMatRef);
42
43			%more params
44			setDiagsParams=[];
45			if iscell(setDiags);
46			setDiagsParams={setDiags{2:end}};
47			setDiags=setDiags{1};
48			end;
49
50			%%%%%%%%%%%%%%%%%%%%%%
51			%load grid and params:
52			%%%%%%%%%%%%%%%%%%%%%%
53
54			gcmfaces_global; global myparms;
55			test1=~isempty(dir([dirMat 'basic_diags_ecco_mygrid.mat']));
56			test2=~isempty(dir([dirMat 'diags_grid_parms.mat']));
57			if ~test1&~test2;
58			error('missing diags_grid_parms.mat')
59			elseif test2;
60			nameGrid='diags_grid_parms.mat';
61			suffDiag='diags_set_';
62			budgetList='diags_select_budget_list.mat';
63			else;
64			nameGrid='basic_diags_ecco_mygrid.mat';
65			suffDiag='basic_diags_ecco_';
66			budgetList='basic_diags_ecco_budget_list.mat';
67			end;
68
69			%here we always reload the grid from dirMat to make sure the same one is used throughout
70			eval(['load ' dirMat nameGrid ';']);
71
72	gforget	1.2	%backward compatibility:
73			if ~isfield(mygrid,'memoryLimit'); mygrid.memoryLimit=0; end;
74			if ~isfield(mygrid,'ioSize'); mygrid.ioSize=0; end;
75
76	gforget	1.1	%in case mygrid.memoryLimit=1, load the stuff that was not saved to diags_grid_parms.mat
77			if mygrid.memoryLimit==1;
78			list0={'hFacS','hFacW'};
79			for iFld=1:length(list0);
80			eval(['mygrid.' list0{iFld} '=rdmds2gcmfaces([mygrid.dirGrid ''' list0{iFld} '*'']);']);
81			end;
82			%
83			mygrid.hFacCsurf=mygrid.hFacC;
84			for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
85			%
86			mskC=mygrid.hFacC; mskC(mskC==0)=NaN; mskC(mskC>0)=1; mygrid.mskC=mskC;
87			mskW=mygrid.hFacW; mskW(mskW==0)=NaN; mskW(mskW>0)=1; mygrid.mskW=mskW;
88			mskS=mygrid.hFacS; mskS(mskS==0)=NaN; mskS(mskS>0)=1; mygrid.mskS=mskS;
89			%
90			gcmfaces_lines_zonal;
91			mygrid.LATS=[mygrid.LATS_MASKS.lat]';
92			[lonPairs,latPairs,names]=line_greatC_TUV_MASKS_v4;
93			gcmfaces_lines_transp(lonPairs,latPairs,names);
94			end;
95
96			%%%%%%%%%%%%%%
97			%define diags:
98			%%%%%%%%%%%%%%
99
100			if strcmp(nameGrid,'diags_grid_parms.mat');
101			%load listDiags and get fileMat (if applies)
102			userStep=1;
103			if ~isempty(which(['diags_set_' setDiags]));
104			eval(['diags_set_' setDiags]);
105			else;
106			diags_set_user;
107			end;
108
109			%reformat listDiags to cell array:
110			jj=strfind(listDiags,' '); jj=[0 jj length(listDiags)+1];
111			for ii=1:length(jj)-1;
112			tmp1=listDiags(jj(ii)+1:jj(ii+1)-1);
113			if ii==1; listDiags2={tmp1}; else; listDiags2{ii}=tmp1; end;
114			end;
115			listDiags=listDiags2; clear listDiags2;
116			end;
117
118			%set fileMat to default name if needed
119			if isempty(who('fileMat')); fileMat=[suffDiag setDiags]; end;
120
121			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
122			%are there results to display?
123			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124
125			test1=isempty(dir([dirMat '/' fileMat '_*.mat']));
126			tmp1=[dirMat '/' fileMat '_*.mat']; tmp2=strfind(tmp1,'_'); tmp1(tmp2)=' ';
127			if test1&addToTex;
128			eval(['load ' dirTex '/myPlots.mat;']);
129			write2tex(fileTex,3,{'[ ',mySection,' ]'},1);
130			write2tex(fileTex,3,{'abort : did not find any'},1);
131			write2tex(fileTex,3,{tmp1},1);
132			write2tex(fileTex,3,{'results files to display'},1);
133			return;
134			elseif test1;
135			fprintf(['\n abort : did not find any \n ' tmp1 '\n results files to display\n\n']);
136			return;
137			end;
138
139
140			%%%%%%%%%%%%%%%%%%%%%%%%%
141			%load pre-computed diags:
142			%%%%%%%%%%%%%%%%%%%%%%%%%
143
144			tic;
145
146			diagsWereLoaded=0;
147
148			%specific load sequence (if any, diagsWereLoaded will be set to 1)
149			if strcmp(nameGrid,'diags_grid_parms.mat');
150			userStep=0;
151			if ~isempty(which(['diags_set_' setDiags]));
152			eval(['diags_set_' setDiags]);
153			else;
154			diags_set_user;
155			end;
156			end;
157
158			%generic load (if no specific one, or )
159			if ~diagsWereLoaded;
160			alldiag=alldiag_load(dirMat,[fileMat '_*.mat']);
161			end;
162
163			toc;
164
165			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
166			%load ref and compute anomalies:
167			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
168
169			if doAnomalies;
170			% scaleAnom=10;
171			scaleAnom=1;
172
173			%store:
174			dirMatBak=dirMat; dirMat=dirMatRef; alldiagBak=alldiag;
175
176			%load reference solution:
177			diagsWereLoaded=0;
178
179			%specific load sequence (if any, diagsWereLoaded will be set to 1)
180			if strcmp(nameGrid,'diags_grid_parms.mat');
181			userStep=0;
182			if ~isempty(which(['diags_set_' setDiags]));
183			eval(['diags_set_' setDiags]);
184			else;
185			diags_set_user;
186			end;
187			end;
188
189			%generic load (if no specific one, or )
190			if ~diagsWereLoaded;
191			alldiag=alldiag_load(dirMat,[fileMat '_*.mat']);
192			end;
193
194			%restore:
195			dirMat=dirMatBak; alldiagRef=alldiag; alldiag=alldiagBak;
196
197			%compute anomalies:
198			for ii=1:length(alldiag.listDiags);
199			tmp0=alldiag.listDiags{ii};
200	gforget	1.2	if ~strcmp(tmp0,'listTimes')&~strcmp(tmp0,'listSteps')
201	gforget	1.1	if isfield(alldiagRef,tmp0);%compute difference
202			tmp1=getfield(alldiag,tmp0)-getfield(alldiagRef,tmp0);
203			alldiag=setfield(alldiag,tmp0,tmp1);
204			else;%cannot compute diff -> set to NaN
205			tmp1=NaN*getfield(alldiag,tmp0);
206			alldiag=setfield(alldiag,tmp0,tmp1);
207			end;
208	gforget	1.2	end;
209	gforget	1.1	end;
210			end;
211
212
213			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
214			%determine time parameters for plots
215			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
216
217			%number of months for runmean
218			if myparms.diagsAreMonthly;
219			myNmean=12; myNmeanTxt=[' -- ' num2str(myNmean*2) ' months low pass filtered'];
220			else;
221			myNmean=0; myNmeanTxt='';
222			end;
223
224			%first and last year of time average
225			myYmean=myparms.yearInAve;
226			if myYmean(1)==myYmean(2); myYmeanTxt=[num2str(myYmean(1)) ' '];
227			else; myYmeanTxt=[num2str(myYmean(1)) '-' num2str(myYmean(2)) ' '];
228			end;
229
230			%to compute time mean/std on the fly
231			tt=[myparms.recInAve(1):min(myparms.recInAve(2),length(alldiag.listTimes))];
232			TT=alldiag.listTimes(tt);
233			nt=length(TT);
234
235			%if only one record, swicth off time series and variance plots
236			multiTimes=1*(myparms.recInAve(2)>myparms.recInAve(1));
237
238			%if llc90 we can plot overturn etc per basin
239			multiBasins=(sum([90 1170]~=mygrid.ioSize)==0);
240
241			%%%%%%%%%%%%%%%
242			%display diags:
243			%%%%%%%%%%%%%%%
244
245			userStep=-1;
246			if ~isempty(which(['diags_set_' setDiags]));
247			eval(['diags_set_' setDiags]);
248			else;
249			diags_set_user;
250			end;
251
252