matlab_class/ecco_v4/cost_altimeter.m

function []=cost_altimeter(dirModel);
%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

%               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/';
dirOut=[dirModel 'mat/']; if isempty(dir(dirOut)); eval(['mkdir ' dirOut ';']); end;
runName=pwd; tmp1=strfind(runName,'/'); runName=runName(tmp1(end)+1:end);

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

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