matlab_class/ecco_v4/cost_altimeter.m

function []=cost_altimeter(dirModel,dirMat);
%object:       compute or plot the various sea level statistics
%              (std model-obs, model, obs, leading to cost function terms)
%inputs:       dirModel is the model run directory
%              dirMat is the directory where diagnozed .mat files will be saved
%                     -> set it to '' to use the default [dirModel 'mat/']
%              doComp is the switch from computation to plot
doComp=1;
if doComp==1;
  doSave=1;
  doPlot=0;
else;
  doSave=0;
  doPlot=1;
end;
doPrint=0;

for doDifObsOrMod=1:3;

if doDifObsOrMod==1; suf='modMobs'; elseif doDifObsOrMod==2; suf='obs'; else; suf='mod'; end;

%%%%%%%%%%%
%load grid:
%%%%%%%%%%%

gcmfaces_global;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;

%%%%%%%%%%%%%%%
%define pathes:
%%%%%%%%%%%%%%%

dirData='/net/altix3700/raid4/gforget/mysetups/ecco_v4/INPUTS_90x50/';
if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat)); eval(['!mkdir ' dirMat ';']); end;
runName=pwd; tmp1=strfind(runName,'/'); runName=runName(tmp1(end)+1:end);

%%%%%%%%%%%%%%%%%
%do computations:
%%%%%%%%%%%%%%%%%

if doComp==0; 
  eval(['load ' dirMat 'cost_altimeter_' suf '.mat myflds;']);
else;

tic;
%mdt cost function term (misfit plot)
dif_mdt=rdmds2gcmfaces([dirModel 'barfiles/mdtdiff_smooth']);
sig_mdt=v4_read_bin([dirData 'sigma_MDT_glob_eccollc.bin'],1,0);
sig_mdt(find(sig_mdt==0))=NaN;
%store:
myflds.dif_mdt=dif_mdt;
myflds.sig_mdt=sig_mdt;

%skip blanks:
listRecs=[1:7:6147];

toc; tic;
%lsc cost function term:
if doDifObsOrMod==1;
  sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
elseif doDifObsOrMod==2;
  sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
else;
  sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
  sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
end;
%mask missing points:
sladiff_smooth(sladiff_smooth==0)=NaN;
%compute rms:
rms_sladiff_smooth=sqrt(nanmean(sladiff_smooth.^2,3));
std_sladiff_smooth=nanstd(sladiff_smooth,0,3);
%get weight:
sig_sladiff_smooth=v4_read_bin([dirData 'sigma_SLA_smooth_eccollc.bin'],1,0);
sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;
%store:
myflds.rms_sladiff_smooth=rms_sladiff_smooth;
myflds.std_sladiff_smooth=std_sladiff_smooth;
myflds.sig_sladiff_smooth=sig_sladiff_smooth;

toc; tic;
%pointwise/point35days cost function term:
if doDifObsOrMod==1;
  sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
elseif doDifObsOrMod==2;
  sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
else;
  sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
  sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
end;
%mask missing points:
sladiff_point35d(sladiff_point35d==0)=NaN;
%compute rms:
rms_sladiff_point35d=sqrt(nanmean(sladiff_point35d.^2,3));
std_sladiff_point35d=nanstd(sladiff_point35d,0,3);
%store:
myflds.rms_sladiff_point35d=rms_sladiff_point35d;
myflds.std_sladiff_point35d=std_sladiff_point35d;

toc; tic;
%pointwise/1day terms:
sum_all=0*mygrid.XC; msk_all=sum_all;
for ii=1:3; 
  if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
  %topex pointwise misfits:
  if doDifObsOrMod==1;
    sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
  elseif doDifObsOrMod==2;
    sladiff_point=cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
  else;
    sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
    sladiff_point=sladiff_point+cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
  end;
  %compute rms:
  msk_tmp=1*(sladiff_point~=0); 
  msk_tmp=sum(msk_tmp,3); sum_tmp=sum(sladiff_point.^2,3);
  sum_all=sum_all+sum_tmp; msk_all=msk_all+msk_tmp; 
  msk_tmp(find(msk_tmp==0))=NaN;
  eval(['rms_' myset '=sqrt(sum_tmp./msk_tmp);']);
  sladiff_point(sladiff_point==0)=NaN;
  eval(['std_' myset '=nanstd(sladiff_point,0,3);']);
end;
%compute overall rms:
msk_all(find(msk_all==0))=NaN;
rms_sladiff_point=sqrt(sum_all./msk_all);
%fill blanks:
warning('off','MATLAB:divideByZero');
msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>66))=NaN;
rms_sladiff_point=diffsmooth2D_extrap_inv(rms_sladiff_point,msk);
warning('on','MATLAB:divideByZero');
%get weight:
sig_sladiff_point=v4_read_bin([dirData 'sigma_SLA_PWS07r2_glob_eccollc.bin'],1,0);
sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
%store:
myflds.rms_tp=rms_tp; myflds.rms_gfo=rms_gfo; myflds.rms_ers=rms_ers; 
myflds.std_tp=std_tp; myflds.std_gfo=std_gfo; myflds.std_ers=std_ers;
myflds.rms_sladiff_point=rms_sladiff_point;
myflds.sig_sladiff_point=sig_sladiff_point;

%compute zonal mean/median:
for ii=1:4;
  switch ii;
    case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
    case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
    case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
    case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
  end;    
  cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
  cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
end;

toc; %write to disk:
if doSave; eval(['save ' dirMat 'cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;

end;%if doComp

if doPlot; 
cc=[-0.4:0.05:-0.25 -0.2:0.03:-0.05 -0.03:0.01:0.03 0.05:0.03:0.2 0.25:0.05:0.4];
figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
cc=[0:0.005:0.02 0.03:0.01:0.05 0.06:0.02:0.1 0.14:0.03:0.2 0.25:0.05:0.4];
figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow; 
end;

if doPlot&doPrint;
  dirFig='../figs/altimeter/'; ff0=gcf-4;
  for ff=1:4; 
    figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig'); 
    eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);   
    eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']); 
  end;
end;

end;%for doDifObsOrMod=1:3;


function [fldOut]=cost_altimeter_read(fileIn,recIn);

nrec=length(recIn);
global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
lrec=siz(1)*siz(2)*4;
myprec='float32';

fldOut=zeros(siz);
fid=fopen([fileIn '.data'],'r','b');
for irec=1:nrec;
status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
end;

fldOut=convert2gcmfaces(fldOut);


1	gforget	1.7	function []=cost_altimeter(dirModel,dirMat);
2			%object: compute or plot the various sea level statistics
3			% (std model-obs, model, obs, leading to cost function terms)
4			%inputs: dirModel is the model run directory
5			% dirMat is the directory where diagnozed .mat files will be saved
6			% -> set it to '' to use the default [dirModel 'mat/']
7			% doComp is the switch from computation to plot
8	gforget	1.4	doComp=1;
9			if doComp==1;
10			doSave=1;
11			doPlot=0;
12			else;
13			doSave=0;
14			doPlot=1;
15			end;
16			doPrint=0;
17
18			for doDifObsOrMod=1:3;
19	gforget	1.1
20	gforget	1.4	if doDifObsOrMod==1; suf='modMobs'; elseif doDifObsOrMod==2; suf='obs'; else; suf='mod'; end;
21	gforget	1.1
22			%%%%%%%%%%%
23			%load grid:
24			%%%%%%%%%%%
25
26	gforget	1.6	gcmfaces_global;
27			if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
28			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
29	gforget	1.1
30			%%%%%%%%%%%%%%%
31			%define pathes:
32			%%%%%%%%%%%%%%%
33
34			dirData='/net/altix3700/raid4/gforget/mysetups/ecco_v4/INPUTS_90x50/';
35	gforget	1.7	if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
36	gforget	1.8	if isempty(dir(dirMat)); eval(['!mkdir ' dirMat ';']); end;
37	gforget	1.1	runName=pwd; tmp1=strfind(runName,'/'); runName=runName(tmp1(end)+1:end);
38
39			%%%%%%%%%%%%%%%%%
40			%do computations:
41			%%%%%%%%%%%%%%%%%
42
43			if doComp==0;
44	gforget	1.7	eval(['load ' dirMat 'cost_altimeter_' suf '.mat myflds;']);
45	gforget	1.1	else;
46
47			tic;
48			%mdt cost function term (misfit plot)
49	gforget	1.4	dif_mdt=rdmds2gcmfaces([dirModel 'barfiles/mdtdiff_smooth']);
50	gforget	1.1	sig_mdt=v4_read_bin([dirData 'sigma_MDT_glob_eccollc.bin'],1,0);
51	gforget	1.5	sig_mdt(find(sig_mdt==0))=NaN;
52	gforget	1.1	%store:
53			myflds.dif_mdt=dif_mdt;
54			myflds.sig_mdt=sig_mdt;
55
56	gforget	1.4	%skip blanks:
57			listRecs=[1:7:6147];
58
59	gforget	1.1	toc; tic;
60			%lsc cost function term:
61	gforget	1.4	if doDifObsOrMod==1;
62			sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
63			elseif doDifObsOrMod==2;
64			sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
65			else;
66			sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
67			sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
68			end;
69			%mask missing points:
70			sladiff_smooth(sladiff_smooth==0)=NaN;
71	gforget	1.1	%compute rms:
72			rms_sladiff_smooth=sqrt(nanmean(sladiff_smooth.^2,3));
73	gforget	1.9	std_sladiff_smooth=nanstd(sladiff_smooth,0,3);
74	gforget	1.1	%get weight:
75			sig_sladiff_smooth=v4_read_bin([dirData 'sigma_SLA_smooth_eccollc.bin'],1,0);
76	gforget	1.5	sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;
77	gforget	1.1	%store:
78			myflds.rms_sladiff_smooth=rms_sladiff_smooth;
79	gforget	1.9	myflds.std_sladiff_smooth=std_sladiff_smooth;
80	gforget	1.1	myflds.sig_sladiff_smooth=sig_sladiff_smooth;
81
82			toc; tic;
83	gforget	1.4	%pointwise/point35days cost function term:
84			if doDifObsOrMod==1;
85			sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
86			elseif doDifObsOrMod==2;
87			sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
88			else;
89			sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
90			sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
91			end;
92			%mask missing points:
93			sladiff_point35d(sladiff_point35d==0)=NaN;
94	gforget	1.1	%compute rms:
95	gforget	1.4	rms_sladiff_point35d=sqrt(nanmean(sladiff_point35d.^2,3));
96	gforget	1.9	std_sladiff_point35d=nanstd(sladiff_point35d,0,3);
97	gforget	1.1	%store:
98	gforget	1.4	myflds.rms_sladiff_point35d=rms_sladiff_point35d;
99	gforget	1.9	myflds.std_sladiff_point35d=std_sladiff_point35d;
100	gforget	1.1
101			toc; tic;
102			%pointwise/1day terms:
103			sum_all=0*mygrid.XC; msk_all=sum_all;
104			for ii=1:3;
105			if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
106			%topex pointwise misfits:
107	gforget	1.4	if doDifObsOrMod==1;
108			sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
109			elseif doDifObsOrMod==2;
110			sladiff_point=cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
111			else;
112			sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
113			sladiff_point=sladiff_point+cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
114			end;
115	gforget	1.1	%compute rms:
116			msk_tmp=1*(sladiff_point~=0);
117			msk_tmp=sum(msk_tmp,3); sum_tmp=sum(sladiff_point.^2,3);
118	gforget	1.5	sum_all=sum_all+sum_tmp; msk_all=msk_all+msk_tmp;
119			msk_tmp(find(msk_tmp==0))=NaN;
120	gforget	1.9	eval(['rms_' myset '=sqrt(sum_tmp./msk_tmp);']);
121			sladiff_point(sladiff_point==0)=NaN;
122			eval(['std_' myset '=nanstd(sladiff_point,0,3);']);
123	gforget	1.1	end;
124			%compute overall rms:
125	gforget	1.5	msk_all(find(msk_all==0))=NaN;
126	gforget	1.1	rms_sladiff_point=sqrt(sum_all./msk_all);
127			%fill blanks:
128	gforget	1.5	warning('off','MATLAB:divideByZero');
129			msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>66))=NaN;
130			rms_sladiff_point=diffsmooth2D_extrap_inv(rms_sladiff_point,msk);
131			warning('on','MATLAB:divideByZero');
132	gforget	1.1	%get weight:
133			sig_sladiff_point=v4_read_bin([dirData 'sigma_SLA_PWS07r2_glob_eccollc.bin'],1,0);
134	gforget	1.5	sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
135	gforget	1.1	%store:
136			myflds.rms_tp=rms_tp; myflds.rms_gfo=rms_gfo; myflds.rms_ers=rms_ers;
137	gforget	1.9	myflds.std_tp=std_tp; myflds.std_gfo=std_gfo; myflds.std_ers=std_ers;
138	gforget	1.1	myflds.rms_sladiff_point=rms_sladiff_point;
139			myflds.sig_sladiff_point=sig_sladiff_point;
140
141			%compute zonal mean/median:
142	gforget	1.4	for ii=1:4;
143	gforget	1.1	switch ii;
144			case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
145			case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
146	gforget	1.4	case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
147			case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
148	gforget	1.1	end;
149	gforget	1.3	cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
150			cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
151	gforget	1.1	end;
152
153			toc; %write to disk:
154	gforget	1.7	if doSave; eval(['save ' dirMat 'cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;
155	gforget	1.1
156			end;%if doComp
157
158			if doPlot;
159	gforget	1.4	cc=[-0.4:0.05:-0.25 -0.2:0.03:-0.05 -0.03:0.01:0.03 0.05:0.03:0.2 0.25:0.05:0.4];
160			figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
161			cc=[0:0.005:0.02 0.03:0.01:0.05 0.06:0.02:0.1 0.14:0.03:0.2 0.25:0.05:0.4];
162			figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
163			figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
164			figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow;
165	gforget	1.1	end;
166
167			if doPlot&doPrint;
168			dirFig='../figs/altimeter/'; ff0=gcf-4;
169			for ff=1:4;
170	gforget	1.4	figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig');
171			eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);
172			eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']);
173	gforget	1.1	end;
174			end;
175
176	gforget	1.4	end;%for doDifObsOrMod=1:3;
177
178
179
180			function [fldOut]=cost_altimeter_read(fileIn,recIn);
181
182			nrec=length(recIn);
183			global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
184			lrec=siz(1)siz(2)4;
185			myprec='float32';
186
187			fldOut=zeros(siz);
188			fid=fopen([fileIn '.data'],'r','b');
189			for irec=1:nrec;
190			status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
191			fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
192			end;
193
194			fldOut=convert2gcmfaces(fldOut);
195
196
197