matlab_class/ecco_v4/cost_sst.m

function []=cost_sst(dirModel,doComp,varargin);
%object:        compute cost function term for sst 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');

%
choiceWeights=2;
choiceData=1;

if choiceWeights==1;
fld_err=read2memory('/net/weddell/raid3/gforget/ecco_v4/input_files/sigma_SST_eccollc.bin',[90 1170]);
fld_err(find(fld_err==0))=NaN;
fld_err(find(fld_err<0.5))=0.5; 
fprintf('I need to fix the error field...\n');
elseif choiceWeights==2;
fld_err=read2memory('/net/weddell/raid3/gforget/ecco_v4/input_files/sigma_half.bin',[90 1170]);
end;

fld_err=convert2gcmfaces(fld_err);
fld_w=fld_err.^-2;

dirData='/net/weddell/raid3/gforget/ecco_v4/input_files/';
fileModel=dir([dirModel 'barfiles/tbar*data']); fileModel=['barfiles/' fileModel.name];

fld_dif=convert2gcmfaces(NaN*ones(90,1170,12*17));
for ycur=1992:2007;
for mcur=1:12;
  if choiceData==1;
    fld_dat=v4_read_bin([dirData 'gael_quick_interp/g_SST_monthly_r2_' num2str(ycur)],mcur,0);
  else;
    fld_dat=v4_read_bin([dirData 'SST_REYNOLDS/REYNOLDS_gael_' num2str(ycur)],mcur,0);
  end;
  fld_dat(find(fld_dat==0))=NaN;
  fld_dat(find(fld_dat<-99))=NaN;

  mm=(ycur-1992)*12+mcur;
  fld_mod=v4_read_bin([dirModel fileModel],mm,1).*mygrid.mskC(:,:,1);
 
  fld_dif(:,:,mm)=fld_mod-fld_dat;
end;
end;

fld_rms=sqrt(nanmean(fld_dif.^2,3));
fld_cost=fld_rms.^2.*fld_w;

eval(['save ' dirModel '/mat/cost_sst.mat fld_err fld_rms fld_cost;']);

else;%display previously computed results

global mygrid;

eval(['load ' dirModel '/mat/cost_sst.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_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});

figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
myCaption={'rms modeled-observed -- sea surface temperature (K)'};
if addToTex; write2tex([dirModel 'tex/myPlots.tex'],2,myCaption,gcf); end;

end;


1	gforget	1.1	function []=cost_sst(dirModel,doComp,varargin);
2			%object: compute cost function term for sst 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			%
17			choiceWeights=2;
18			choiceData=1;
19
20			if choiceWeights==1;
21			fld_err=read2memory('/net/weddell/raid3/gforget/ecco_v4/input_files/sigma_SST_eccollc.bin',[90 1170]);
22			fld_err(find(fld_err==0))=NaN;
23			fld_err(find(fld_err<0.5))=0.5;
24			fprintf('I need to fix the error field...\n');
25			elseif choiceWeights==2;
26			fld_err=read2memory('/net/weddell/raid3/gforget/ecco_v4/input_files/sigma_half.bin',[90 1170]);
27			end;
28
29			fld_err=convert2gcmfaces(fld_err);
30			fld_w=fld_err.^-2;
31
32			dirData='/net/weddell/raid3/gforget/ecco_v4/input_files/';
33			fileModel=dir([dirModel 'barfiles/tbar*data']); fileModel=['barfiles/' fileModel.name];
34
35			fld_dif=convert2gcmfaces(NaNones(90,1170,1217));
36			for ycur=1992:2007;
37			for mcur=1:12;
38			if choiceData==1;
39			fld_dat=v4_read_bin([dirData 'gael_quick_interp/g_SST_monthly_r2_' num2str(ycur)],mcur,0);
40			else;
41			fld_dat=v4_read_bin([dirData 'SST_REYNOLDS/REYNOLDS_gael_' num2str(ycur)],mcur,0);
42			end;
43			fld_dat(find(fld_dat==0))=NaN;
44			fld_dat(find(fld_dat<-99))=NaN;
45
46			mm=(ycur-1992)*12+mcur;
47			fld_mod=v4_read_bin([dirModel fileModel],mm,1).*mygrid.mskC(:,:,1);
48
49			fld_dif(:,:,mm)=fld_mod-fld_dat;
50			end;
51			end;
52
53			fld_rms=sqrt(nanmean(fld_dif.^2,3));
54			fld_cost=fld_rms.^2.*fld_w;
55
56			eval(['save ' dirModel '/mat/cost_sst.mat fld_err fld_rms fld_cost;']);
57
58			else;%display previously computed results
59
60			global mygrid;
61
62			eval(['load ' dirModel '/mat/cost_sst.mat;']);
63
64			%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
65			%figure; m_map_gcmfaces(fld_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
66
67			figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
68			myCaption={'rms modeled-observed -- sea surface temperature (K)'};
69			if addToTex; write2tex([dirModel 'tex/myPlots.tex'],2,myCaption,gcf); end;
70
71			end;
72
73