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));
%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 ' 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			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			%get weight:
74			sig_sladiff_smooth=v4_read_bin([dirData 'sigma_SLA_smooth_eccollc.bin'],1,0);
75	gforget	1.5	sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;
76	gforget	1.1	%store:
77			myflds.rms_sladiff_smooth=rms_sladiff_smooth;
78			myflds.sig_sladiff_smooth=sig_sladiff_smooth;
79
80			toc; tic;
81	gforget	1.4	%pointwise/point35days cost function term:
82			if doDifObsOrMod==1;
83			sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
84			elseif doDifObsOrMod==2;
85			sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
86			else;
87			sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
88			sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
89			end;
90			%mask missing points:
91			sladiff_point35d(sladiff_point35d==0)=NaN;
92	gforget	1.1	%compute rms:
93	gforget	1.4	rms_sladiff_point35d=sqrt(nanmean(sladiff_point35d.^2,3));
94	gforget	1.1	%store:
95	gforget	1.4	myflds.rms_sladiff_point35d=rms_sladiff_point35d;
96	gforget	1.1
97			toc; tic;
98			%pointwise/1day terms:
99			sum_all=0*mygrid.XC; msk_all=sum_all;
100			for ii=1:3;
101			if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
102			%topex pointwise misfits:
103	gforget	1.4	if doDifObsOrMod==1;
104			sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
105			elseif doDifObsOrMod==2;
106			sladiff_point=cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
107			else;
108			sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
109			sladiff_point=sladiff_point+cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
110			end;
111	gforget	1.1	%compute rms:
112			msk_tmp=1*(sladiff_point~=0);
113			msk_tmp=sum(msk_tmp,3); sum_tmp=sum(sladiff_point.^2,3);
114	gforget	1.5	sum_all=sum_all+sum_tmp; msk_all=msk_all+msk_tmp;
115			msk_tmp(find(msk_tmp==0))=NaN;
116	gforget	1.1	eval(['rms_' myset '=sum_tmp./msk_tmp;']);
117			end;
118			%compute overall rms:
119	gforget	1.5	msk_all(find(msk_all==0))=NaN;
120	gforget	1.1	rms_sladiff_point=sqrt(sum_all./msk_all);
121			%fill blanks:
122	gforget	1.5	warning('off','MATLAB:divideByZero');
123			msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>66))=NaN;
124			rms_sladiff_point=diffsmooth2D_extrap_inv(rms_sladiff_point,msk);
125			warning('on','MATLAB:divideByZero');
126	gforget	1.1	%get weight:
127			sig_sladiff_point=v4_read_bin([dirData 'sigma_SLA_PWS07r2_glob_eccollc.bin'],1,0);
128	gforget	1.5	sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
129	gforget	1.1	%store:
130			myflds.rms_tp=rms_tp; myflds.rms_gfo=rms_gfo; myflds.rms_ers=rms_ers;
131			myflds.rms_sladiff_point=rms_sladiff_point;
132			myflds.sig_sladiff_point=sig_sladiff_point;
133
134			%compute zonal mean/median:
135	gforget	1.4	for ii=1:4;
136	gforget	1.1	switch ii;
137			case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
138			case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
139	gforget	1.4	case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
140			case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
141	gforget	1.1	end;
142	gforget	1.3	cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
143			cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
144	gforget	1.1	end;
145
146			toc; %write to disk:
147	gforget	1.7	if doSave; eval(['save ' dirMat 'cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;
148	gforget	1.1
149			end;%if doComp
150
151			if doPlot;
152	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];
153			figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
154			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];
155			figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
156			figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
157			figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow;
158	gforget	1.1	end;
159
160			if doPlot&doPrint;
161			dirFig='../figs/altimeter/'; ff0=gcf-4;
162			for ff=1:4;
163	gforget	1.4	figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig');
164			eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);
165			eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']);
166	gforget	1.1	end;
167			end;
168
169	gforget	1.4	end;%for doDifObsOrMod=1:3;
170
171
172
173			function [fldOut]=cost_altimeter_read(fileIn,recIn);
174
175			nrec=length(recIn);
176			global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
177			lrec=siz(1)siz(2)4;
178			myprec='float32';
179
180			fldOut=zeros(siz);
181			fid=fopen([fileIn '.data'],'r','b');
182			for irec=1:nrec;
183			status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
184			fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
185			end;
186
187			fldOut=convert2gcmfaces(fldOut);
188
189
190