matlab_class/ecco_v4/cost_sst.m

function []=cost_sst(dirModel,dirMat,doComp,dirTex,nameTex);
%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)
%               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;

%grid, params and inputs

gcmfaces_global; global myparms;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
if isfield(myparms,'yearFirst'); yearFirst=myparms.yearFirst; yearLast=myparms.yearLast;
else; yearFirst=1992; yearLast=2011;
end;

dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';
subdirReynolds='input_sst_reynolds_v4/'; nameReynolds='reynolds_oiv2_r1';
subdirRemss='input_sst_remss_v4/'; nameRemss='tmi_amsre_oisst_r1';
subdirL2p='input_sst_amsre_daily/'; nameL2p='tmi_amsre_l2p_r1';

if 0;%old versions
nameReynolds='v4_reynolds_monthly';
nameRemss='v4_TMI_AMSRE_sst_monav';
nameL2p='amsre_r2';
end;

dirErr='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_all_from_pleiades/';
fld_err=read2memory([dirErr '/sigma_half.bin'],[90 1170]);
fld_err=convert2gcmfaces(fld_err);
fld_w=fld_err.^-2;

fileModel=dir([dirModel 'barfiles/tbar*data']); fileModel=['barfiles/' fileModel.name];

%computational loop
mod_m_rey=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
mod_m_remss=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
%
zm_mod=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_rey=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_mod_m_rey=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_remss=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_mod_m_remss=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
%
for ycur=yearFirst:yearLast;
fprintf(['starting ' num2str(ycur) '\n']);
tic;
for mcur=1:12;
  %load Reynolds SST
  file0=[dirData subdirReynolds nameReynolds '_' num2str(ycur)];
  if ~isempty(dir(file0)); fld_rey=v4_read_bin(file0,mcur,0); else; fld_rey=NaN*mygrid.mskC(:,:,1); end;
  fld_rey(find(fld_rey==0))=NaN;
  fld_rey(find(fld_rey<-99))=NaN;

  %load Remss SST
  file0=[dirData subdirRemss nameRemss '_' num2str(ycur)];
  if ~isempty(dir(file0)); fld_remss=v4_read_bin(file0,mcur,0); else; fld_remss=NaN*mygrid.mskC(:,:,1); end;
  fld_remss(find(fld_remss==0))=NaN;
  fld_remss(find(fld_remss<-99))=NaN;

  %load model SST
  mm=(ycur-yearFirst)*12+mcur;
  fld_mod=v4_read_bin([dirModel fileModel],mm,1).*mygrid.mskC(:,:,1);

  %store misfit maps 
  mod_m_rey(:,:,mm)=fld_mod-fld_rey;
  mod_m_remss(:,:,mm)=fld_mod-fld_remss;

  %comp zonal means
  zm_mod(:,mm)=calc_zonmean_T(fld_mod.*mygrid.mskC(:,:,1));

  msk=mygrid.mskC(:,:,1); msk(fld_rey==0)=NaN;
  zm_rey(:,mm)=calc_zonmean_T(fld_rey.*msk);
  zm_mod_m_rey(:,mm)=calc_zonmean_T(fld_mod-fld_rey.*msk);

  msk=mygrid.mskC(:,:,1); msk(fld_remss==0)=NaN;
  zm_remss(:,mm)=calc_zonmean_T(fld_remss.*msk);
  zm_mod_m_remss(:,mm)=calc_zonmean_T(fld_mod-fld_remss.*msk);
end;
toc;
end;

%compute rms maps
rms_to_rey=sqrt(nanmean(mod_m_rey.^2,3));
rms_to_remss=sqrt(nanmean(mod_m_remss.^2,3));

eval(['save ' dirMat '/cost_sst.mat fld_err rms_* zm_*;']);

else;%display previously computed results

global mygrid;

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

if isempty(who('fld_rms')); fld_rms=rms_to_rey; end;

%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.6	function []=cost_sst(dirModel,dirMat,doComp,dirTex,nameTex);
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.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	gforget	1.7	%grid, params and inputs
28
29			gcmfaces_global; global myparms;
30	gforget	1.3	if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
31			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
32	gforget	1.7	if isfield(myparms,'yearFirst'); yearFirst=myparms.yearFirst; yearLast=myparms.yearLast;
33			else; yearFirst=1992; yearLast=2011;
34			end;
35	gforget	1.1
36	gforget	1.7	dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';
37			subdirReynolds='input_sst_reynolds_v4/'; nameReynolds='reynolds_oiv2_r1';
38			subdirRemss='input_sst_remss_v4/'; nameRemss='tmi_amsre_oisst_r1';
39			subdirL2p='input_sst_amsre_daily/'; nameL2p='tmi_amsre_l2p_r1';
40
41			if 0;%old versions
42			nameReynolds='v4_reynolds_monthly';
43			nameRemss='v4_TMI_AMSRE_sst_monav';
44			nameL2p='amsre_r2';
45	gforget	1.1	end;
46
47	gforget	1.7	dirErr='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_all_from_pleiades/';
48			fld_err=read2memory([dirErr '/sigma_half.bin'],[90 1170]);
49	gforget	1.1	fld_err=convert2gcmfaces(fld_err);
50			fld_w=fld_err.^-2;
51
52			fileModel=dir([dirModel 'barfiles/tbar*data']); fileModel=['barfiles/' fileModel.name];
53
54	gforget	1.7	%computational loop
55			mod_m_rey=repmat(NaNmygrid.mskC(:,:,1),[1 1 12(yearLast-yearFirst+1)]);
56			mod_m_remss=repmat(NaNmygrid.mskC(:,:,1),[1 1 12(yearLast-yearFirst+1)]);
57			%
58			zm_mod=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
59			zm_rey=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
60			zm_mod_m_rey=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
61			zm_remss=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
62			zm_mod_m_remss=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
63			%
64			for ycur=yearFirst:yearLast;
65			fprintf(['starting ' num2str(ycur) '\n']);
66			tic;
67	gforget	1.1	for mcur=1:12;
68	gforget	1.7	%load Reynolds SST
69			file0=[dirData subdirReynolds nameReynolds '_' num2str(ycur)];
70			if ~isempty(dir(file0)); fld_rey=v4_read_bin(file0,mcur,0); else; fld_rey=NaN*mygrid.mskC(:,:,1); end;
71			fld_rey(find(fld_rey==0))=NaN;
72			fld_rey(find(fld_rey<-99))=NaN;
73
74			%load Remss SST
75			file0=[dirData subdirRemss nameRemss '_' num2str(ycur)];
76			if ~isempty(dir(file0)); fld_remss=v4_read_bin(file0,mcur,0); else; fld_remss=NaN*mygrid.mskC(:,:,1); end;
77			fld_remss(find(fld_remss==0))=NaN;
78			fld_remss(find(fld_remss<-99))=NaN;
79	gforget	1.1
80	gforget	1.7	%load model SST
81			mm=(ycur-yearFirst)*12+mcur;
82	gforget	1.1	fld_mod=v4_read_bin([dirModel fileModel],mm,1).*mygrid.mskC(:,:,1);
83	gforget	1.7
84			%store misfit maps
85			mod_m_rey(:,:,mm)=fld_mod-fld_rey;
86			mod_m_remss(:,:,mm)=fld_mod-fld_remss;
87
88			%comp zonal means
89			zm_mod(:,mm)=calc_zonmean_T(fld_mod.*mygrid.mskC(:,:,1));
90
91			msk=mygrid.mskC(:,:,1); msk(fld_rey==0)=NaN;
92			zm_rey(:,mm)=calc_zonmean_T(fld_rey.*msk);
93			zm_mod_m_rey(:,mm)=calc_zonmean_T(fld_mod-fld_rey.*msk);
94
95			msk=mygrid.mskC(:,:,1); msk(fld_remss==0)=NaN;
96			zm_remss(:,mm)=calc_zonmean_T(fld_remss.*msk);
97			zm_mod_m_remss(:,mm)=calc_zonmean_T(fld_mod-fld_remss.*msk);
98	gforget	1.1	end;
99	gforget	1.7	toc;
100	gforget	1.1	end;
101
102	gforget	1.7	%compute rms maps
103			rms_to_rey=sqrt(nanmean(mod_m_rey.^2,3));
104			rms_to_remss=sqrt(nanmean(mod_m_remss.^2,3));
105	gforget	1.1
106	gforget	1.7	eval(['save ' dirMat '/cost_sst.mat fld_err rms_* zm_*;']);
107	gforget	1.1
108			else;%display previously computed results
109
110			global mygrid;
111
112	gforget	1.4	eval(['load ' dirMat '/cost_sst.mat;']);
113	gforget	1.1
114	gforget	1.7	if isempty(who('fld_rms')); fld_rms=rms_to_rey; end;
115
116	gforget	1.1	%figure; m_map_gcmfaces(fld_cost,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]});
117			%figure; m_map_gcmfaces(fld_err,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
118
119			figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
120	gforget	1.2	myCaption={'modeled-observed rms -- sea surface temperature (K)'};
121	gforget	1.4	if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
122	gforget	1.1
123			end;
124
125