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;

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

if ~isempty(dir([dirModel 'sigma_MDT_glob_eccollc.bin']));
  dirSigma=dirModel;
else;
  dirSigma='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_err/';
end;

if isempty(dir([dirModel 'barfiles']));
  dirEcco=dirModel;
else;
  dirEcco=[dirModel 'barfiles' filesep];
end;

%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_r4.err','slaerr_gridscale_r4.err'}; maxLatObs=90

if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat));     mkdir([dirMat]); end;
if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
runName=pwd; tmp1=strfind(runName,filesep); runName=runName(tmp1(end)+1:end);

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

for doDifObsOrMod=1:3;

    if doDifObsOrMod==1; suf='modMobs'; elseif doDifObsOrMod==2; suf='obs'; else; suf='mod'; end;
    
    if doComp==0;
        eval(['load ' dirMat 'cost/cost_altimeter_' suf '.mat myflds;']);
    else;
        
        tic;
        %mdt cost function term (misfit plot)
        dif_mdt=rdmds2gcmfaces([dirEcco 'mdtdiff_smooth']);
        sig_mdt=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([dirEcco 'sladiff_smooth*data']);
        nrec=tmp1.bytes/90/1170/4;
        ttShift=17;
        %skip 1992
        listRecs=[1+365-ttShift:17+365*18-ttShift];
        nRecs=length(listRecs); TT=1993+[0:nRecs-1]/365.25;
        
        toc; tic;
        %pointwise/1day terms:
        sladiff_point=repmat(0*mygrid.RAC,[1 1 nRecs]); count_point=sladiff_point;
        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_tmp=cost_altimeter_read([dirEcco 'sladiff_' myset '_raw'],ttShift+listRecs);
            elseif doDifObsOrMod==2;
                sladiff_tmp=cost_altimeter_read([dirEcco 'slaobs_' myset '_raw'],ttShift+listRecs);
            else;
                sladiff_tmp=cost_altimeter_read([dirEcco 'sladiff_' myset '_raw'],ttShift+listRecs);
                sladiff_tmp=sladiff_tmp+cost_altimeter_read([dirEcco 'slaobs_' myset '_raw'],ttShift+listRecs);
            end;
            %add to overall data set:
            sladiff_point=sladiff_point+sladiff_tmp; count_point=count_point+(sladiff_tmp~=0);
            %finalize data set:
            sladiff_tmp(sladiff_tmp==0)=NaN;
            %compute rms:
            count_tmp=sum(~isnan(sladiff_tmp),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
            eval(['myflds.rms_' myset '=sqrt(nansum(sladiff_tmp.^2,3)./count_tmp);']);
            eval(['myflds.std_' myset '=nanstd(sladiff_tmp,0,3).*msk_tmp;']);
        end;
        %finalize overall data set:
        count_point(count_point==0)=NaN;
        sladiff_point=sladiff_point./count_point;
        %compute overall rms,std:
        count_tmp=sum(~isnan(sladiff_point),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
        myflds.rms_sladiff_point=sqrt(nansum(sladiff_point.^2,3)./count_tmp);
        myflds.std_sladiff_point=nanstd(sladiff_point,0,3).*msk_tmp;
        %
        msk_point=1*(count_point>0);
        %fill blanks:
        warning('off','MATLAB:divideByZero');
        msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>maxLatObs))=NaN;
        myflds.xtrp_rms_sladiff_point=diffsmooth2D_extrap_inv(myflds.rms_sladiff_point,msk);
        myflds.xtrp_std_sladiff_point=diffsmooth2D_extrap_inv(myflds.std_sladiff_point,msk);
        warning('on','MATLAB:divideByZero');
        %get weight:
        sig_sladiff_point=read_bin([dirSigma nameSigma{3}],1,0);
        sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
        myflds.sig_sladiff_point=sig_sladiff_point;
        
        %computational mask : only points that were actually observed, and in 35d average
        msk_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
        msk_point35d=1*(msk_point35d~=0);
        tmp1=sum(msk_point35d==0&msk_point~=0);
        tmp2=sum(msk_point35d~=0&msk_point~=0);
        fprintf('after masking : %d omitted, %d retained \n',tmp1,tmp2);
        msk_point=msk_point.*msk_point35d;
        msk_point(msk_point==0)=NaN;
        %plotting mask : regions of less than 100 observations are omitted
        msk_100pts=1*(nansum(msk_point,3)>=100);
        msk_100pts(msk_100pts==0)=NaN;
        %store
        myflds.msk_100pts=msk_100pts;
        
        toc; tic;
        %lsc cost function term:
        if doDifObsOrMod==1;
            sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
        elseif doDifObsOrMod==2;
            sladiff_smooth=cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
        else;
            sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
            sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
        end;
        %mask missing points:
        sladiff_smooth(sladiff_smooth==0)=NaN;
        sladiff_smooth=sladiff_smooth.*msk_point;
        count_tmp=sum(~isnan(sladiff_smooth),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
        %compute rms:
        rms_sladiff_smooth=msk_tmp.*sqrt(nanmean(sladiff_smooth.^2,3));
        std_sladiff_smooth=msk_tmp.*nanstd(sladiff_smooth,0,3);
        %get weight:
        sig_sladiff_smooth=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([dirEcco 'sladiff_raw'],listRecs);
        elseif doDifObsOrMod==2;
            sladiff_point35d=cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
        else;
            sladiff_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
            sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
        end;
        %mask missing points:
        sladiff_point35d(sladiff_point35d==0)=NaN;
        sladiff_point35d=sladiff_point35d.*msk_point;
        count_tmp=sum(~isnan(sladiff_point35d),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
        %compute rms:
        rms_sladiff_point35d=msk_tmp.*sqrt(nanmean(sladiff_point35d.^2,3));
        std_sladiff_point35d=msk_tmp.*nanstd(sladiff_point35d,0,3);
        %store:
        myflds.rms_sladiff_point35d=rms_sladiff_point35d;
        myflds.std_sladiff_point35d=std_sladiff_point35d;
        
        if 0;
            %difference between scales
            rms_sladiff_35dMsmooth=msk_tmp.*sqrt(nanmean((sladiff_point35d-sladiff_smooth).^2,3));
            std_sladiff_35dMsmooth=msk_tmp.*nanstd((sladiff_point35d-sladiff_smooth),0,3);
            %store:
            myflds.rms_sladiff_35dMsmooth=rms_sladiff_35dMsmooth;
            myflds.std_sladiff_35dMsmooth=std_sladiff_35dMsmooth;

            %difference between scales
            rms_sladiff_pointMpoint35d=msk_tmp.*sqrt(nanmean((sladiff_point-sladiff_point35d).^2,3));
            std_sladiff_pointMpoint35d=msk_tmp.*nanstd((sladiff_point-sladiff_point35d),0,3);
            %store:
            myflds.rms_sladiff_pointMpoint35d=rms_sladiff_pointMpoint35d;
            myflds.std_sladiff_pointMpoint35d=std_sladiff_pointMpoint35d;
        end;
        
        %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/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.13
9	gforget	1.4	doComp=1;
10			if doComp==1;
11	gforget	1.11	doSave=1;
12			doPlot=0;
13	gforget	1.4	else;
14	gforget	1.11	doSave=0;
15			doPlot=1;
16	gforget	1.4	end;
17			doPrint=0;
18
19	gforget	1.14	%%%%%%%%%%%
20			%load grid:
21			%%%%%%%%%%%
22
23			gcmfaces_global;
24			if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
25			if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
26
27			%%%%%%%%%%%%%%%
28			%define pathes:
29			%%%%%%%%%%%%%%%
30
31			if ~isempty(dir([dirModel 'sigma_MDT_glob_eccollc.bin']));
32			dirSigma=dirModel;
33			else;
34			dirSigma='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_err/';
35			end;
36
37			if isempty(dir([dirModel 'barfiles']));
38			dirEcco=dirModel;
39			else;
40			dirEcco=[dirModel 'barfiles' filesep];
41			end;
42
43			%nameSigma={'sigma_MDT_glob_eccollc.bin','sigma_SLA_smooth_eccollc.bin','sigma_SLA_PWS07r2_glob_eccollc.bin'}; maxLatObs=66
44			nameSigma={'sigma_MDT_glob_eccollc.bin','slaerr_largescale_r4.err','slaerr_gridscale_r4.err'}; maxLatObs=90
45
46			if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
47			if isempty(dir(dirMat)); mkdir([dirMat]); end;
48			if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
49			runName=pwd; tmp1=strfind(runName,filesep); runName=runName(tmp1(end)+1:end);
50
51			%%%%%%%%%%%%%%%%%
52			%do computations:
53			%%%%%%%%%%%%%%%%%
54
55	gforget	1.4	for doDifObsOrMod=1:3;
56	gforget	1.14
57	gforget	1.11	if doDifObsOrMod==1; suf='modMobs'; elseif doDifObsOrMod==2; suf='obs'; else; suf='mod'; end;
58
59			if doComp==0;
60			eval(['load ' dirMat 'cost/cost_altimeter_' suf '.mat myflds;']);
61			else;
62
63			tic;
64			%mdt cost function term (misfit plot)
65	gforget	1.14	dif_mdt=rdmds2gcmfaces([dirEcco 'mdtdiff_smooth']);
66	gforget	1.12	sig_mdt=read_bin([dirSigma nameSigma{1}],1,0);
67	gforget	1.11	sig_mdt(find(sig_mdt==0))=NaN;
68			%store:
69			myflds.dif_mdt=dif_mdt;
70			myflds.sig_mdt=sig_mdt;
71
72			%skip blanks:
73	gforget	1.14	tmp1=dir([dirEcco 'sladiff_smooth*data']);
74	gforget	1.11	nrec=tmp1.bytes/90/1170/4;
75			ttShift=17;
76			%skip 1992
77			listRecs=[1+365-ttShift:17+365*18-ttShift];
78			nRecs=length(listRecs); TT=1993+[0:nRecs-1]/365.25;
79
80			toc; tic;
81			%pointwise/1day terms:
82			sladiff_point=repmat(0*mygrid.RAC,[1 1 nRecs]); count_point=sladiff_point;
83			for ii=1:3;
84			if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
85			%topex pointwise misfits:
86			if doDifObsOrMod==1;
87	gforget	1.14	sladiff_tmp=cost_altimeter_read([dirEcco 'sladiff_' myset '_raw'],ttShift+listRecs);
88	gforget	1.11	elseif doDifObsOrMod==2;
89	gforget	1.14	sladiff_tmp=cost_altimeter_read([dirEcco 'slaobs_' myset '_raw'],ttShift+listRecs);
90	gforget	1.11	else;
91	gforget	1.14	sladiff_tmp=cost_altimeter_read([dirEcco 'sladiff_' myset '_raw'],ttShift+listRecs);
92			sladiff_tmp=sladiff_tmp+cost_altimeter_read([dirEcco 'slaobs_' myset '_raw'],ttShift+listRecs);
93	gforget	1.11	end;
94			%add to overall data set:
95			sladiff_point=sladiff_point+sladiff_tmp; count_point=count_point+(sladiff_tmp~=0);
96			%finalize data set:
97			sladiff_tmp(sladiff_tmp==0)=NaN;
98			%compute rms:
99			count_tmp=sum(~isnan(sladiff_tmp),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
100			eval(['myflds.rms_' myset '=sqrt(nansum(sladiff_tmp.^2,3)./count_tmp);']);
101			eval(['myflds.std_' myset '=nanstd(sladiff_tmp,0,3).*msk_tmp;']);
102			end;
103			%finalize overall data set:
104			count_point(count_point==0)=NaN;
105			sladiff_point=sladiff_point./count_point;
106			%compute overall rms,std:
107			count_tmp=sum(~isnan(sladiff_point),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
108			myflds.rms_sladiff_point=sqrt(nansum(sladiff_point.^2,3)./count_tmp);
109			myflds.std_sladiff_point=nanstd(sladiff_point,0,3).*msk_tmp;
110			%
111			msk_point=1*(count_point>0);
112			%fill blanks:
113			warning('off','MATLAB:divideByZero');
114			msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>maxLatObs))=NaN;
115			myflds.xtrp_rms_sladiff_point=diffsmooth2D_extrap_inv(myflds.rms_sladiff_point,msk);
116			myflds.xtrp_std_sladiff_point=diffsmooth2D_extrap_inv(myflds.std_sladiff_point,msk);
117			warning('on','MATLAB:divideByZero');
118			%get weight:
119	gforget	1.12	sig_sladiff_point=read_bin([dirSigma nameSigma{3}],1,0);
120	gforget	1.11	sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
121			myflds.sig_sladiff_point=sig_sladiff_point;
122
123			%computational mask : only points that were actually observed, and in 35d average
124	gforget	1.14	msk_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
125	gforget	1.11	msk_point35d=1*(msk_point35d~=0);
126			tmp1=sum(msk_point35d==0&msk_point~=0);
127			tmp2=sum(msk_point35d~=0&msk_point~=0);
128			fprintf('after masking : %d omitted, %d retained \n',tmp1,tmp2);
129			msk_point=msk_point.*msk_point35d;
130			msk_point(msk_point==0)=NaN;
131			%plotting mask : regions of less than 100 observations are omitted
132			msk_100pts=1*(nansum(msk_point,3)>=100);
133			msk_100pts(msk_100pts==0)=NaN;
134			%store
135			myflds.msk_100pts=msk_100pts;
136
137			toc; tic;
138			%lsc cost function term:
139			if doDifObsOrMod==1;
140	gforget	1.14	sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
141	gforget	1.11	elseif doDifObsOrMod==2;
142	gforget	1.14	sladiff_smooth=cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
143	gforget	1.11	else;
144	gforget	1.14	sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
145			sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
146	gforget	1.11	end;
147			%mask missing points:
148			sladiff_smooth(sladiff_smooth==0)=NaN;
149			sladiff_smooth=sladiff_smooth.*msk_point;
150			count_tmp=sum(~isnan(sladiff_smooth),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
151			%compute rms:
152			rms_sladiff_smooth=msk_tmp.*sqrt(nanmean(sladiff_smooth.^2,3));
153			std_sladiff_smooth=msk_tmp.*nanstd(sladiff_smooth,0,3);
154			%get weight:
155	gforget	1.12	sig_sladiff_smooth=read_bin([dirSigma nameSigma{2}],1,0);
156	gforget	1.11	sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;
157			%store:
158			myflds.rms_sladiff_smooth=rms_sladiff_smooth;
159			myflds.std_sladiff_smooth=std_sladiff_smooth;
160			myflds.sig_sladiff_smooth=sig_sladiff_smooth;
161
162			toc; tic;
163			%pointwise/point35days cost function term:
164			if doDifObsOrMod==1;
165	gforget	1.14	sladiff_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
166	gforget	1.11	elseif doDifObsOrMod==2;
167	gforget	1.14	sladiff_point35d=cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
168	gforget	1.11	else;
169	gforget	1.14	sladiff_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
170			sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
171	gforget	1.11	end;
172			%mask missing points:
173			sladiff_point35d(sladiff_point35d==0)=NaN;
174			sladiff_point35d=sladiff_point35d.*msk_point;
175			count_tmp=sum(~isnan(sladiff_point35d),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
176			%compute rms:
177			rms_sladiff_point35d=msk_tmp.*sqrt(nanmean(sladiff_point35d.^2,3));
178			std_sladiff_point35d=msk_tmp.*nanstd(sladiff_point35d,0,3);
179			%store:
180			myflds.rms_sladiff_point35d=rms_sladiff_point35d;
181			myflds.std_sladiff_point35d=std_sladiff_point35d;
182
183			if 0;
184			%difference between scales
185			rms_sladiff_35dMsmooth=msk_tmp.*sqrt(nanmean((sladiff_point35d-sladiff_smooth).^2,3));
186			std_sladiff_35dMsmooth=msk_tmp.*nanstd((sladiff_point35d-sladiff_smooth),0,3);
187			%store:
188			myflds.rms_sladiff_35dMsmooth=rms_sladiff_35dMsmooth;
189			myflds.std_sladiff_35dMsmooth=std_sladiff_35dMsmooth;
190
191			%difference between scales
192			rms_sladiff_pointMpoint35d=msk_tmp.*sqrt(nanmean((sladiff_point-sladiff_point35d).^2,3));
193			std_sladiff_pointMpoint35d=msk_tmp.*nanstd((sladiff_point-sladiff_point35d),0,3);
194			%store:
195			myflds.rms_sladiff_pointMpoint35d=rms_sladiff_pointMpoint35d;
196			myflds.std_sladiff_pointMpoint35d=std_sladiff_pointMpoint35d;
197			end;
198
199			%compute zonal mean/median:
200			for ii=1:4;
201			switch ii;
202			case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
203			case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
204			case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
205			case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
206			end;
207			cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
208			cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
209			end;
210
211			toc; %write to disk:
212			if doSave; eval(['save ' dirMat 'cost/cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;
213
214			end;%if doComp
215
216			if doPlot;
217			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];
218			figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
219			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];
220			figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
221			figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
222			figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow;
223			end;
224
225			if doPlot&doPrint;
226			dirFig='../figs/altimeter/'; ff0=gcf-4;
227			for ff=1:4;
228			figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig');
229			eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);
230			eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']);
231			end;
232			end;
233
234	gforget	1.4	end;%for doDifObsOrMod=1:3;
235
236
237
238			function [fldOut]=cost_altimeter_read(fileIn,recIn);
239
240			nrec=length(recIn);
241			global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
242			lrec=siz(1)siz(2)4;
243			myprec='float32';
244
245			fldOut=zeros(siz);
246			fid=fopen([fileIn '.data'],'r','b');
247			for irec=1:nrec;
248	gforget	1.11	status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
249			fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
250	gforget	1.4	end;
251
252			fldOut=convert2gcmfaces(fldOut);
253
254
255