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:
%%%%%%%%%%%%%%%

dirSigma='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_all_from_pleiades/';
%nameSigma={'sigma_MDT_glob_eccollc.bin','sigma_SLA_smooth_eccollc.bin','sigma_SLA_PWS07r2_glob_eccollc.bin'}; maxLatObs=66
nameSigma={'sigma_MDT_glob_eccollc.bin','slaerr_largescale_r1.err','slaerr_gridscale_r1.err'}; maxLatObs=90

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([dirSigma nameSigma{1}],1,0);
sig_mdt(find(sig_mdt==0))=NaN;
%store:
myflds.dif_mdt=dif_mdt;
myflds.sig_mdt=sig_mdt;

%skip blanks:
tmp1=dir([dirModel 'barfiles/sladiff_smooth*data']);
nrec=tmp1.bytes/90/1170/4;
listRecs=[1:7:nrec];

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