matlab_class/gcmfaces_diags/diags_display.m

function [alldiag]=diags_display(dirMat,setDiags,dirTex,nameTex);
% [alldiag]=DIAGS_DISPLAY(dirMat,setDiags,dirTex,nameTex)
%    displays a set of diagnostics (setDiags) from the results
%    stored (in [dirMat 'diags_set_' setDiags]). 
%
%    Further outputs the results to tex ([dirTex nameTex '.tex'])
%    if dirTex and nameTex are specified. This functionality 
%    is normally operated via diags_driver_tex.m. If dirMat is 
%    specified as {dirMat,dirMatRef} then anomalies will be plotted. 
%
%    setDiags is the choice of diagnostics set such as
%                       'A') trasnports
%                       'B') air-sea fluxes
%                       'C') state variables
%                       'D') global and hemispheric budgets
%                       'MLD') mixed layer depths
%
%    In some cases more specific options can be specified, e.g. to 
%    display only subsurface budgets one may set setDiags={{'D',11}}
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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=1; fileTex='';
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;

%reload myparms from dirMat (and mygrid if included the mat file)
eval(['load ' dirMat nameGrid ';']);

%reload mygrid if needed
if isfield(myparms,'dirGrid'); diags_grid(myparms.dirGrid,0); end;

%zonal mean and sections needed for transport computations
if ~isfield(mygrid,'mygrid.LATS_MASKS');
    gcmfaces_lines_zonal;
end;
if ~isfield(mygrid,'LINES_MASKS');
    [lonPairs,latPairs,names]=gcmfaces_lines_pairs;
    gcmfaces_lines_transp(lonPairs,latPairs,names);
end;

%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]=gcmfaces_lines_pairs;
        gcmfaces_lines_transp(lonPairs,latPairs,names);
end;

%%%%%%%%%%%%%%
%define diags:
%%%%%%%%%%%%%%

%test for results organization
if isdir([dirMat 'diags_set_' setDiags]); 
    dirMat=[dirMat 'diags_set_' setDiags '/']; 
end;
if isdir([dirMatRef 'diags_set_' setDiags]); 
    dirMatRef=[dirMatRef 'diags_set_' setDiags '/']; 
end;

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