matlab_class/ecco_v4/cost_bp.m

function []=cost_bp(dirModel,doComp,varargin);
%object:        compute cost function term for grace data
%inputs:        dimodel is the model directory
%               doComp is a switch (1->compute; 0->display)
%optional:      addToTex

if nargin>2; addToTex=varargin{1}; else; addToTex=0; end;

if doComp;

%load grid
dirGrid='GRID/';
global mygrid;
gcmfaces_path; grid_load(dirGrid,5,'compact');

%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/weddell/raid3/gforget/ecco_v4/input_files/grace_tmp/GRACE_CSR_v4.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 ' dirModel '/mat/cost_bp.mat fld_err fld_dif fld_mod fld_dat fld_cost;']);

else;%display previously computed results

global mygrid;

eval(['load ' dirModel '/mat/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={'rms modeled-observed -- bottom pressure (cm)'};
if addToTex; write2tex([dirModel 'tex/myPlots.tex'],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([dirModel 'tex/myPlots.tex'],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([dirModel 'tex/myPlots.tex'],2,myCaption,gcf); end;

end;

1	gforget	1.1	function []=cost_bp(dirModel,doComp,varargin);
2			%object: compute cost function term for grace data
3			%inputs: dimodel is the model directory
4			% doComp is a switch (1->compute; 0->display)
5			%optional: addToTex
6
7			if nargin>2; addToTex=varargin{1}; else; addToTex=0; end;
8
9			if doComp;
10
11			%load grid
12			dirGrid='GRID/';
13			global mygrid;
14			gcmfaces_path; grid_load(dirGrid,5,'compact');
15
16			%read model cost output
17			fld_dat=rdmds2gcmfaces([dirModel 'barfiles/bpdatanom_smooth']);
18			fld_dif=rdmds2gcmfaces([dirModel 'barfiles/bpdifanom_smooth']);
19
20			%mask:
21			fld_msk=rdmds2gcmfaces([dirModel 'barfiles/bpdatanom_raw']);
22			fld_msk(find(fld_msk~=0))=1;
23			fld_msk(find(fld_msk==0))=NaN;
24
25			fld_dif=fld_dif.*fld_msk;
26			fld_dat=fld_dat.*fld_msk;
27			fld_mod=fld_dat+fld_dif;%compute model values
28
29			%read uncertainty fields
30			fld_err=read2memory('/net/weddell/raid3/gforget/ecco_v4/input_files/grace_tmp/GRACE_CSR_v4.err',[90 1170]);
31			fld_err=convert2gcmfaces(fld_err);
32			fld_err(find(fld_err==0))=NaN;
33
34			%compute weight
35			fld_w=fld_err.^-2;
36
37			%compute cost
38			fld_cost=mk3d(fld_w,fld_dif).*(fld_dif.^2);
39			fld_cost=nanmean(fld_cost,3);
40
41			fld_dif=nanstd(fld_dif,[],3);
42			fld_mod=nanstd(fld_mod,[],3);
43			fld_dat=nanstd(fld_dat,[],3);
44
45			eval(['save ' dirModel '/mat/cost_bp.mat fld_err fld_dif fld_mod fld_dat fld_cost;']);
46
47			else;%display previously computed results
48
49			global mygrid;
50
51			eval(['load ' dirModel '/mat/cost_bp.mat;']);
52
53			%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
54			figure; m_map_gcmfaces(fld_dif,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
55			myCaption={'rms modeled-observed -- bottom pressure (cm)'};
56			if addToTex; write2tex([dirModel 'tex/myPlots.tex'],2,myCaption,gcf); end;
57
58			figure; m_map_gcmfaces(fld_mod,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
59			myCaption={'rms modeled -- bottom pressure (cm)'};
60			if addToTex; write2tex([dirModel 'tex/myPlots.tex'],2,myCaption,gcf); end;
61
62			figure; m_map_gcmfaces(fld_dat,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
63			myCaption={'rms observed -- bottom pressure (cm)'};
64			if addToTex; write2tex([dirModel 'tex/myPlots.tex'],2,myCaption,gcf); end;
65
66			end;
67