matlab_class/ecco_v4/cost_bp.m

function []=cost_bp(dirModel,dirMat,doComp,dirTex,nameTex);
%object:        compute cost function term for grace data
%inputs:        dimodel 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/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat)); eval(['!mkdir ' dirMat ';']); end;

%determine if and where to create tex and figures files
dirMat=[dirMat '/'];
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 nameTex '.tex'];
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;

%read model cost output
fld_dat=rdmds2gcmfaces([dirModel 'barfiles/bpdatanom_smooth']);
fld_dif=rdmds2gcmfaces([dirModel 'barfiles/bpdifanom_smooth']);

%mask:
fld_msk=rdmds2gcmfaces([dirModel 'barfiles/bpdatanom_raw']);
fld_msk(find(fld_msk~=0))=1; 
fld_msk(find(fld_msk==0))=NaN;

fld_dif=fld_dif.*fld_msk;
fld_dat=fld_dat.*fld_msk;
fld_mod=fld_dat+fld_dif;%compute model values

%read uncertainty fields
fld_err=read2memory('/net/nares/raid11/ecco-shared/ecco-version-4/input/input_grace/GRACE_CSR_withland_err',[90 1170]);
fld_err=convert2gcmfaces(fld_err);
fld_err(find(fld_err==0))=NaN;

%compute weight
fld_w=fld_err.^-2;

%compute cost
fld_cost=mk3d(fld_w,fld_dif).*(fld_dif.^2); 
fld_cost=nanmean(fld_cost,3);

fld_dif=nanstd(fld_dif,[],3);
fld_mod=nanstd(fld_mod,[],3);
fld_dat=nanstd(fld_dat,[],3);

eval(['save ' dirMat '/cost_bp.mat fld_err fld_dif fld_mod fld_dat fld_cost;']);

else;%display previously computed results

global mygrid;

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

%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
figure; m_map_gcmfaces(fld_dif,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
myCaption={'modeled-observed rms -- bottom pressure (cm)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure; m_map_gcmfaces(fld_mod,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
myCaption={'rms modeled -- bottom pressure (cm)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure; m_map_gcmfaces(fld_dat,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
myCaption={'rms observed -- bottom pressure (cm)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
myCaption={'Cost function'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

end;

1	gforget	1.6	function []=cost_bp(dirModel,dirMat,doComp,dirTex,nameTex);
2	gforget	1.1	%object: compute cost function term for grace data
3			%inputs: dimodel is the model directory
4	gforget	1.4	% 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.4	%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.6	% nameTex is the tex file name (default : 'myPlots')
10	gforget	1.1
11	gforget	1.4	if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
12	gforget	1.5	if isempty(dir(dirMat)); eval(['!mkdir ' dirMat ';']); end;
13	gforget	1.4
14			%determine if and where to create tex and figures files
15			dirMat=[dirMat '/'];
16			if isempty(who('dirTex'));
17			addToTex=0;
18			else;
19			if ~ischar(dirTex); error('mis-specified dirTex'); end;
20	gforget	1.6	addToTex=1;
21			if isempty(who('nameTex')); nameTex='myPlots'; end;
22			fileTex=[dirTex nameTex '.tex'];
23	gforget	1.4	end;
24	gforget	1.1
25			if doComp;
26
27			%load grid
28	gforget	1.3	gcmfaces_global;
29			if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
30			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
31	gforget	1.1
32			%read model cost output
33			fld_dat=rdmds2gcmfaces([dirModel 'barfiles/bpdatanom_smooth']);
34			fld_dif=rdmds2gcmfaces([dirModel 'barfiles/bpdifanom_smooth']);
35
36			%mask:
37			fld_msk=rdmds2gcmfaces([dirModel 'barfiles/bpdatanom_raw']);
38			fld_msk(find(fld_msk~=0))=1;
39			fld_msk(find(fld_msk==0))=NaN;
40
41			fld_dif=fld_dif.*fld_msk;
42			fld_dat=fld_dat.*fld_msk;
43			fld_mod=fld_dat+fld_dif;%compute model values
44
45			%read uncertainty fields
46	gforget	1.7	fld_err=read2memory('/net/nares/raid11/ecco-shared/ecco-version-4/input/input_grace/GRACE_CSR_withland_err',[90 1170]);
47	gforget	1.1	fld_err=convert2gcmfaces(fld_err);
48			fld_err(find(fld_err==0))=NaN;
49
50			%compute weight
51			fld_w=fld_err.^-2;
52
53			%compute cost
54			fld_cost=mk3d(fld_w,fld_dif).*(fld_dif.^2);
55			fld_cost=nanmean(fld_cost,3);
56
57			fld_dif=nanstd(fld_dif,[],3);
58			fld_mod=nanstd(fld_mod,[],3);
59			fld_dat=nanstd(fld_dat,[],3);
60
61	gforget	1.4	eval(['save ' dirMat '/cost_bp.mat fld_err fld_dif fld_mod fld_dat fld_cost;']);
62	gforget	1.1
63			else;%display previously computed results
64
65			global mygrid;
66
67	gforget	1.4	eval(['load ' dirMat '/cost_bp.mat;']);
68	gforget	1.1
69			%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
70			figure; m_map_gcmfaces(fld_dif,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
71	gforget	1.2	myCaption={'modeled-observed rms -- bottom pressure (cm)'};
72	gforget	1.4	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
73	gforget	1.1
74			figure; m_map_gcmfaces(fld_mod,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
75			myCaption={'rms modeled -- bottom pressure (cm)'};
76	gforget	1.4	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
77	gforget	1.1
78			figure; m_map_gcmfaces(fld_dat,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
79			myCaption={'rms observed -- bottom pressure (cm)'};
80	gforget	1.4	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
81	gforget	1.1
82	gforget	1.7	figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
83			myCaption={'Cost function'};
84			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
85
86	gforget	1.1	end;
87