matlab_class/ecco_v4/cost_sst.m

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

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; fileTex=[dirTex 'myPlots.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;

%
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 ' dirMat '/cost_sst.mat fld_err fld_rms fld_cost;']);

else;%display previously computed results

global mygrid;

eval(['load ' dirMat '/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={'modeled-observed rms -- sea surface temperature (K)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

end;


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