matlab_class/ecco_v4/cost_xx.m

function []=cost_xx(dirModel,dirMat,doComp,dirTex,nameTex);
%object:        compute cost function term for atmospheric controls
%inputs:        dirModel is the model directory
%               dirMat is the directory where diagnozed .mat files will be saved
%                     -> set it to '' to use the default [dirModel 'mat/']
%               doComp is a switch (1->compute; 0->display)
%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')

if isempty(dirMat); dirMat=[dirModel 'mat' filesep]; else; dirMat=[dirMat filesep]; end;
if isempty(dir(dirMat));     mkdir([dirMat]); end;

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

for ii=1:6;
switch ii;
case 1; xxName='atemp'; nameSigma='cap_sigma_tmp2m_degC_eccollc.bin'; cc=2; uni='K';
case 2; xxName='aqh'; nameSigma='cap_sigma_spfh2m_eccollc.bin'; cc=2; uni='g/kg';
case 3; xxName='tauu'; nameSigma='cap_sigma_ustr_eccollc.bin'; cc=0.04; uni='N/m2';
case 4; xxName='tauv'; nameSigma='cap_sigma_vstr_eccollc.bin'; cc=0.04; uni='N/m2';
case 5; xxName='lwdown'; nameSigma='cap_sigma_dlw_eccollc.bin'; cc=20; uni='W/m2';
case 6; xxName='swdown'; nameSigma='cap_sigma_dsw_eccollc.bin'; cc=40; uni='W/m2';
end;

if doComp;

%load grid
gcmfaces_global;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;

if ~isempty(dir([dirModel nameSigma])); 
  dirSig=dirModel;
else;
  warning(['could not find ' nameSigma ' -> abort']);
  return;
end;

if ~isempty(dir([dirModel 'xx_' xxName '.effective.*data']));
  dirXX=dirModel;
else;
  dirXX=[dirModel 'ADXXfiles' filesep];
end;

%read model cost output
tmp1=dir([dirXX 'xx_' xxName '.effective.*data']);
tmp2=size(convert2gcmfaces(mygrid.XC)); 
fld_xx=read2memory([dirXX tmp1.name],[tmp2 tmp1.bytes/tmp2(1)/tmp2(2)/4]);
fld_xx=convert2gcmfaces(fld_xx);
%does not work when adjoint overwrites xx*effective*meta     fld_xx=rdmds2gcmfaces([dirXX 'xx_' xxName '.effective.*']);

fld_sig=read_bin([dirSig nameSigma],1,0);
if strcmp(xxName,'aqh'); fld_xx=fld_xx*1000; fld_sig=fld_sig*1000; end;

%compute xx stats
fld_rms=sqrt(mean(fld_xx.^2,3));
fld_mean=mean(fld_xx,3);
fld_std=std(fld_xx,[],3);

%mask
fld_rms=fld_rms.*mygrid.mskC(:,:,1); 
fld_sig=fld_sig.*mygrid.mskC(:,:,1);
fld_mean=fld_mean.*mygrid.mskC(:,:,1);
fld_std=fld_std.*mygrid.mskC(:,:,1);

clear fld_xx;

if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
eval(['save ' dirMat '/cost/cost_xx_' xxName '.mat fld_* cc uni;']);

else;%display previously computed results

global mygrid;

if ~isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;

eval(['load ' dirMat '/cost_xx_' xxName '.mat;']);

figure; 
m_map_gcmfaces(fld_sig,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
myCaption={['prior uncertainty -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure; 
m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
myCaption={['rms adjustment -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure;
m_map_gcmfaces(fld_std,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
myCaption={['std adjustment -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure;
m_map_gcmfaces(fld_mean,0,{'myCaxis',[-0.5:0.05:0.5]*cc});
myCaption={['mean adjustment -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

end;

end;%for ii=1:6;

1	gforget	1.5	function []=cost_xx(dirModel,dirMat,doComp,dirTex,nameTex);
2	gforget	1.1	%object: compute cost function term for atmospheric controls
3			%inputs: dirModel is the model directory
4	gforget	1.3	% dirMat is the directory where diagnozed .mat files will be saved
5			% -> set it to '' to use the default [dirModel 'mat/']
6	gforget	1.1	% doComp is a switch (1->compute; 0->display)
7	gforget	1.3	%optional: dirTex is the directory where tex and figures files are created
8			% (if not specified then display all results to screen instead)
9	gforget	1.5	% nameTex is the tex file name (default : 'myPlots')
10	gforget	1.1
11	gforget	1.11	if isempty(dirMat); dirMat=[dirModel 'mat' filesep]; else; dirMat=[dirMat filesep]; end;
12	gforget	1.9	if isempty(dir(dirMat)); mkdir([dirMat]); end;
13	gforget	1.3
14			%determine if and where to create tex and figures files
15	gforget	1.11	dirMat=[dirMat filesep];
16	gforget	1.3	if isempty(who('dirTex'));
17			addToTex=0;
18			else;
19			if ~ischar(dirTex); error('mis-specified dirTex'); end;
20	gforget	1.5	addToTex=1;
21			if isempty(who('nameTex')); nameTex='myPlots'; end;
22	gforget	1.11	fileTex=[dirTex filesep nameTex '.tex'];
23	gforget	1.3	end;
24	gforget	1.1
25			for ii=1:6;
26			switch ii;
27	gforget	1.12	case 1; xxName='atemp'; nameSigma='cap_sigma_tmp2m_degC_eccollc.bin'; cc=2; uni='K';
28			case 2; xxName='aqh'; nameSigma='cap_sigma_spfh2m_eccollc.bin'; cc=2; uni='g/kg';
29			case 3; xxName='tauu'; nameSigma='cap_sigma_ustr_eccollc.bin'; cc=0.04; uni='N/m2';
30			case 4; xxName='tauv'; nameSigma='cap_sigma_vstr_eccollc.bin'; cc=0.04; uni='N/m2';
31			case 5; xxName='lwdown'; nameSigma='cap_sigma_dlw_eccollc.bin'; cc=20; uni='W/m2';
32			case 6; xxName='swdown'; nameSigma='cap_sigma_dsw_eccollc.bin'; cc=40; uni='W/m2';
33	gforget	1.1	end;
34
35			if doComp;
36
37			%load grid
38	gforget	1.2	gcmfaces_global;
39			if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
40			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
41
42	gforget	1.12	if ~isempty(dir([dirModel nameSigma]));
43	gforget	1.10	dirSig=dirModel;
44			else;
45	gforget	1.13	warning(['could not find ' nameSigma ' -> abort']);
46			return;
47	gforget	1.10	end;
48
49			if ~isempty(dir([dirModel 'xx_' xxName '.effective.*data']));
50			dirXX=dirModel;
51			else;
52			dirXX=[dirModel 'ADXXfiles' filesep];
53			end;
54	gforget	1.1
55			%read model cost output
56	gforget	1.10	tmp1=dir([dirXX 'xx_' xxName '.effective.*data']);
57	gforget	1.4	tmp2=size(convert2gcmfaces(mygrid.XC));
58	gforget	1.10	fld_xx=read2memory([dirXX tmp1.name],[tmp2 tmp1.bytes/tmp2(1)/tmp2(2)/4]);
59	gforget	1.4	fld_xx=convert2gcmfaces(fld_xx);
60	gforget	1.10	%does not work when adjoint overwrites xxeffectivemeta fld_xx=rdmds2gcmfaces([dirXX 'xx_' xxName '.effective.*']);
61	gforget	1.4
62	gforget	1.12	fld_sig=read_bin([dirSig nameSigma],1,0);
63	gforget	1.1	if strcmp(xxName,'aqh'); fld_xx=fld_xx1000; fld_sig=fld_sig1000; end;
64
65			%compute xx stats
66			fld_rms=sqrt(mean(fld_xx.^2,3));
67			fld_mean=mean(fld_xx,3);
68			fld_std=std(fld_xx,[],3);
69
70			%mask
71			fld_rms=fld_rms.*mygrid.mskC(:,:,1);
72			fld_sig=fld_sig.*mygrid.mskC(:,:,1);
73			fld_mean=fld_mean.*mygrid.mskC(:,:,1);
74			fld_std=fld_std.*mygrid.mskC(:,:,1);
75
76			clear fld_xx;
77
78	gforget	1.7	if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
79			eval(['save ' dirMat '/cost/cost_xx_' xxName '.mat fld_* cc uni;']);
80	gforget	1.1
81			else;%display previously computed results
82
83			global mygrid;
84
85	gforget	1.7	if ~isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
86
87	gforget	1.3	eval(['load ' dirMat '/cost_xx_' xxName '.mat;']);
88	gforget	1.1
89			figure;
90			m_map_gcmfaces(fld_sig,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
91			myCaption={['prior uncertainty -- ' xxName ' (' uni ')']};
92	gforget	1.3	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
93	gforget	1.1
94			figure;
95			m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
96			myCaption={['rms adjustment -- ' xxName ' (' uni ')']};
97	gforget	1.3	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
98	gforget	1.1
99			figure;
100			m_map_gcmfaces(fld_std,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
101			myCaption={['std adjustment -- ' xxName ' (' uni ')']};
102	gforget	1.3	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
103	gforget	1.1
104			figure;
105			m_map_gcmfaces(fld_mean,0,{'myCaxis',[-0.5:0.05:0.5]*cc});
106			myCaption={['mean adjustment -- ' xxName ' (' uni ')']};
107	gforget	1.3	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
108	gforget	1.1
109			end;
110
111			end;%for ii=1:6;
112