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' filesep]; else; dirMat=[dirMat filesep]; end;
if isempty(dir(dirMat));     mkdir([dirMat]); end;

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

%search for nctiles files
useNctiles=0;
dirNctiles=[dirModel 'nctiles_remotesensing/sst/'];
if ~isempty(dir(dirNctiles)); useNctiles=1; end;

nameReynolds='reynolds_oiv2_r1';
nameRemss='tmi_amsre_oisst_r1';
nameL2p='tmi_amsre_l2p_r1';

dirData=dirModel; subdirReynolds=''; subdirRemss=''; subdirL2p='';
if ~isempty(dir([dirModel 'inputfiles/' nameReynolds '*']));
  dirData=[dirModel 'inputfiles/']; subdirReynolds=''; subdirRemss=''; subdirL2p='';
end;

nameSigma='sigma_surf_0p5.bin';
dirSigma=dirModel;
if ~isempty(dir([dirModel 'inputfiles/' nameSigma]));
  dirSigma=[dirModel 'inputfiles/'];
end;

fld_err=NaN*mygrid.RAC;
file0=[dirSigma nameSigma];
if ~isempty(dir(file0)); fld_err=read_bin(file0,1,0); end;
if useNctiles; fld_err=read_nctiles([dirNctiles 'sst'],'sst_sigma'); end;
fld_w=fld_err.^-2;

nameTbar='m_sst_mon';
if ~isempty(dir([dirModel nameTbar '*']));
  dirTbar=dirModel;
else;
  dirTbar=[dirModel 'barfiles' filesep];
end;
fileModel='unknown';
file0=[dirTbar nameTbar '*data']; 
if ~isempty(dir(file0)); 
  fileModel=dir(file0);
  fileModel=fileModel.name;
end;

if ~useNctiles;
%computational loop
sst_mod=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
sst_rey=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
sst_remss=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
%
for ycur=yearFirst:yearLast;
fprintf(['starting ' num2str(ycur) '\n']);
tic;
for mcur=1:12;

  mm=(ycur-yearFirst)*12+mcur;

  %load Reynolds SST
  fld_rey=NaN*mygrid.RAC;
  file0=[dirData subdirReynolds nameReynolds '_' num2str(ycur)];
  if ~isempty(dir(file0)); fld_rey=read_bin(file0,mcur,0); end;
  fld_rey(find(fld_rey==0))=NaN;
  fld_rey(find(fld_rey<-99))=NaN;
  sst_rey(:,:,mm)=fld_rey;

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

  %load model SST
  fld_mod=NaN*mygrid.RAC;
  file0=[dirTbar fileModel];
  if ~isempty(dir(file0)); fld_mod=read_bin(file0,mm,1); end;
  fld_mod=fld_mod.*mygrid.mskC(:,:,1);
  sst_mod(:,:,mm)=fld_mod;
end;
toc;
end;
end;%if ~useNctiles;

if useNctiles; 
 sst_rey=read_nctiles([dirNctiles 'sst'],'sst_REYNOLDS');
 sst_mod=read_nctiles([dirNctiles 'sst'],'sst_ECCOv4r2');
 sst_remss=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
end;

%store misfit maps
mod_m_rey=sst_mod-sst_rey;
mod_m_remss=sst_mod-sst_remss;

%comp zonal means
tic;
zm_mod=calc_zonmean_T(sst_mod);
zm_rey=calc_zonmean_T(sst_rey);
zm_remss=calc_zonmean_T(sst_remss);
zm_mod_m_rey=calc_zonmean_T(mod_m_rey);
zm_mod_m_remss=calc_zonmean_T(mod_m_remss);
toc;

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

if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
eval(['save ' dirMat '/cost/cost_sst.mat fld_err rms_* zm_*;']);

else;%display previously computed results

global mygrid;

if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;

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

if ~isempty(who('fld_rms')); 
  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;
else;

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

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

  ny=size(zm_rey,2)/12;
  [xx,yy]=meshgrid(1992+([1:ny*12]-0.5)/12,mygrid.LATS);

  figureL;
  obs=zm_rey; obsCycle=sst_cycle(obs); 
  mod=zm_rey+zm_mod_m_rey; modCycle=sst_cycle(mod);
  mis=zm_mod_m_rey; misCycle=sst_cycle(mis);
  subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar;
  set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude'); 
  title('Reynolds sst anomaly');
  subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar;
  set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
  title('ECCO sst anomaly');
  subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar;
  set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-Reynolds sst misfit');
  myCaption={'ECCO and Reynolds zonal mean sst anomalies (K)'};
  if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

  figureL;
  obs=zm_remss; obsCycle=sst_cycle(obs);
  mod=zm_remss+zm_mod_m_remss; modCycle=sst_cycle(mod);
  mis=zm_mod_m_remss; misCycle=sst_cycle(mis);
  subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar;
  set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
  title('REMSS sst anomaly');
  subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar;
  set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
  title('ECCO sst anomaly');
  subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar;
  set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-REMSS sst misfit');
  myCaption={'ECCO and REMSS zonal mean sst anomalies (K)'};
  if addToTex; write2tex(fileTex,2,myCaption,gcf); end;

end;

end;

function [zmCycle]=sst_cycle(zmIn);

ny=size(zmIn,2)/12;
zmCycle=NaN*zeros(179,12);
for mm=1:12;
zmCycle(:,mm)=nanmean(zmIn(:,mm:12:ny*12),2);
end;
zmCycle=repmat(zmCycle,[1 ny]);

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.14	if isempty(dirMat); dirMat=[dirModel 'mat' filesep]; else; dirMat=[dirMat filesep]; end;
12	gforget	1.12	if isempty(dir(dirMat)); mkdir([dirMat]); end;
13	gforget	1.4
14			%determine if and where to create tex and figures files
15	gforget	1.14	dirMat=[dirMat filesep];
16	gforget	1.4	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	gforget	1.14	fileTex=[dirTex filesep 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.17	%search for nctiles files
37			useNctiles=0;
38			dirNctiles=[dirModel 'nctiles_remotesensing/sst/'];
39			if ~isempty(dir(dirNctiles)); useNctiles=1; end;
40
41	gforget	1.13	nameReynolds='reynolds_oiv2_r1';
42			nameRemss='tmi_amsre_oisst_r1';
43			nameL2p='tmi_amsre_l2p_r1';
44	gforget	1.17
45			dirData=dirModel; subdirReynolds=''; subdirRemss=''; subdirL2p='';
46			if ~isempty(dir([dirModel 'inputfiles/' nameReynolds '*']));
47	gforget	1.15	dirData=[dirModel 'inputfiles/']; subdirReynolds=''; subdirRemss=''; subdirL2p='';
48	gforget	1.1	end;
49
50	gforget	1.16	nameSigma='sigma_surf_0p5.bin';
51	gforget	1.17	dirSigma=dirModel;
52			if ~isempty(dir([dirModel 'inputfiles/' nameSigma]));
53	gforget	1.15	dirSigma=[dirModel 'inputfiles/'];
54	gforget	1.13	end;
55	gforget	1.15
56	gforget	1.17	fld_err=NaN*mygrid.RAC;
57			file0=[dirSigma nameSigma];
58			if ~isempty(dir(file0)); fld_err=read_bin(file0,1,0); end;
59			if useNctiles; fld_err=read_nctiles([dirNctiles 'sst'],'sst_sigma'); end;
60	gforget	1.1	fld_w=fld_err.^-2;
61
62	gforget	1.17	nameTbar='m_sst_mon';
63	gforget	1.13	if ~isempty(dir([dirModel nameTbar '*']));
64			dirTbar=dirModel;
65			else;
66			dirTbar=[dirModel 'barfiles' filesep];
67			end;
68	gforget	1.17	fileModel='unknown';
69			file0=[dirTbar nameTbar '*data'];
70			if ~isempty(dir(file0));
71			fileModel=dir(file0);
72			fileModel=fileModel.name;
73			end;
74	gforget	1.1
75	gforget	1.17	if ~useNctiles;
76	gforget	1.7	%computational loop
77	gforget	1.17	sst_mod=repmat(NaNmygrid.mskC(:,:,1),[1 1 12(yearLast-yearFirst+1)]);
78			sst_rey=repmat(NaNmygrid.mskC(:,:,1),[1 1 12(yearLast-yearFirst+1)]);
79			sst_remss=repmat(NaNmygrid.mskC(:,:,1),[1 1 12(yearLast-yearFirst+1)]);
80	gforget	1.7	%
81			for ycur=yearFirst:yearLast;
82			fprintf(['starting ' num2str(ycur) '\n']);
83			tic;
84	gforget	1.1	for mcur=1:12;
85	gforget	1.17
86			mm=(ycur-yearFirst)*12+mcur;
87
88	gforget	1.7	%load Reynolds SST
89	gforget	1.17	fld_rey=NaN*mygrid.RAC;
90	gforget	1.7	file0=[dirData subdirReynolds nameReynolds '_' num2str(ycur)];
91	gforget	1.17	if ~isempty(dir(file0)); fld_rey=read_bin(file0,mcur,0); end;
92	gforget	1.7	fld_rey(find(fld_rey==0))=NaN;
93			fld_rey(find(fld_rey<-99))=NaN;
94	gforget	1.17	sst_rey(:,:,mm)=fld_rey;
95	gforget	1.7
96			%load Remss SST
97	gforget	1.17	fld_remss=NaN*mygrid.RAC;
98	gforget	1.7	file0=[dirData subdirRemss nameRemss '_' num2str(ycur)];
99	gforget	1.17	if ~isempty(dir(file0)); fld_remss=read_bin(file0,mcur,0); end;
100	gforget	1.7	fld_remss(find(fld_remss==0))=NaN;
101			fld_remss(find(fld_remss<-99))=NaN;
102	gforget	1.17	sst_remss(:,:,mm)=fld_remss;
103	gforget	1.1
104	gforget	1.7	%load model SST
105	gforget	1.17	fld_mod=NaN*mygrid.RAC;
106			file0=[dirTbar fileModel];
107			if ~isempty(dir(file0)); fld_mod=read_bin(file0,mm,1); end;
108			fld_mod=fld_mod.*mygrid.mskC(:,:,1);
109			sst_mod(:,:,mm)=fld_mod;
110			end;
111			toc;
112			end;
113			end;%if ~useNctiles;
114	gforget	1.7
115	gforget	1.17	if useNctiles;
116			sst_rey=read_nctiles([dirNctiles 'sst'],'sst_REYNOLDS');
117			sst_mod=read_nctiles([dirNctiles 'sst'],'sst_ECCOv4r2');
118			sst_remss=repmat(NaNmygrid.mskC(:,:,1),[1 1 12(yearLast-yearFirst+1)]);
119	gforget	1.1	end;
120	gforget	1.17
121			%store misfit maps
122			mod_m_rey=sst_mod-sst_rey;
123			mod_m_remss=sst_mod-sst_remss;
124
125			%comp zonal means
126			tic;
127			zm_mod=calc_zonmean_T(sst_mod);
128			zm_rey=calc_zonmean_T(sst_rey);
129			zm_remss=calc_zonmean_T(sst_remss);
130			zm_mod_m_rey=calc_zonmean_T(mod_m_rey);
131			zm_mod_m_remss=calc_zonmean_T(mod_m_remss);
132	gforget	1.7	toc;
133	gforget	1.1
134	gforget	1.7	%compute rms maps
135			rms_to_rey=sqrt(nanmean(mod_m_rey.^2,3));
136			rms_to_remss=sqrt(nanmean(mod_m_remss.^2,3));
137	gforget	1.1
138	gforget	1.10	if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
139			eval(['save ' dirMat '/cost/cost_sst.mat fld_err rms_* zm_*;']);
140	gforget	1.1
141			else;%display previously computed results
142
143			global mygrid;
144
145	gforget	1.9	if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
146
147	gforget	1.4	eval(['load ' dirMat '/cost_sst.mat;']);
148	gforget	1.1
149	gforget	1.8	if ~isempty(who('fld_rms'));
150			figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
151			myCaption={'modeled-observed rms -- sea surface temperature (K)'};
152			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
153			else;
154
155			figure; m_map_gcmfaces(rms_to_rey,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
156			myCaption={'modeled-Reynolds rms -- sea surface temperature (K)'};
157			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
158
159			figure; m_map_gcmfaces(rms_to_remss,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
160			myCaption={'modeled-REMSS rms -- sea surface temperature (K)'};
161			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
162
163			ny=size(zm_rey,2)/12;
164			[xx,yy]=meshgrid(1992+([1:ny*12]-0.5)/12,mygrid.LATS);
165
166			figureL;
167			obs=zm_rey; obsCycle=sst_cycle(obs);
168			mod=zm_rey+zm_mod_m_rey; modCycle=sst_cycle(mod);
169			mis=zm_mod_m_rey; misCycle=sst_cycle(mis);
170			subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar;
171			set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
172			title('Reynolds sst anomaly');
173			subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar;
174			set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
175			title('ECCO sst anomaly');
176			subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar;
177			set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-Reynolds sst misfit');
178			myCaption={'ECCO and Reynolds zonal mean sst anomalies (K)'};
179			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
180
181			figureL;
182			obs=zm_remss; obsCycle=sst_cycle(obs);
183			mod=zm_remss+zm_mod_m_remss; modCycle=sst_cycle(mod);
184			mis=zm_mod_m_remss; misCycle=sst_cycle(mis);
185			subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar;
186			set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
187			title('REMSS sst anomaly');
188			subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar;
189			set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
190			title('ECCO sst anomaly');
191			subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar;
192			set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-REMSS sst misfit');
193			myCaption={'ECCO and REMSS zonal mean sst anomalies (K)'};
194			if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
195	gforget	1.7
196	gforget	1.8	end;
197
198			end;
199	gforget	1.1
200	gforget	1.8	function [zmCycle]=sst_cycle(zmIn);
201	gforget	1.1
202	gforget	1.8	ny=size(zmIn,2)/12;
203			zmCycle=NaN*zeros(179,12);
204			for mm=1:12;
205			zmCycle(:,mm)=nanmean(zmIn(:,mm:12:ny*12),2);
206	gforget	1.1	end;
207	gforget	1.8	zmCycle=repmat(zmCycle,[1 ny]);
208	gforget	1.1