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

%search for nctiles files
useNctiles=0;
dirNctiles=[dirModel 'nctiles_remotesensing/sealevel/'];
if ~isempty(dir(dirNctiles)); useNctiles=1; end;

nameSigma={'sigma_MDT_glob_eccollc.bin','slaerr_largescale_r5.err','slaerr_gridscale_r5.err'}; maxLatObs=90

dirSigma=dirModel;
if ~isempty(dir([dirModel 'inputfiles/' nameSigma{1}]));
  dirSigma=[dirModel 'inputfiles/'];
end;

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

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);

%%%%%%%%%%%%%%%%%
%global means:
%%%%%%%%%%%%%%%%%

tmpName=dir([dirEcco 'm_eta_day*data']);
avisoName=which('MSL_aviso.bin');
if length(tmpName)==1&~isempty(avisoName);
  m_eta_day=rdmds2gcmfaces([dirEcco tmpName.name(1:end-5)]);
  avisoName=which('MSL_aviso.bin');
  etaglo_aviso=read2memory(avisoName)';
  %global mean
  nt=size(m_eta_day{1},3); etaglo=NaN*zeros(1,nt);
  area=mygrid.mskC(:,:,1).*mygrid.RAC; areatot=nansum(area);
  for tt=1:nt;
    etaglo(tt)=nansum(area.*m_eta_day(:,:,tt))/areatot;
  end;
  %monthly mean
  etaglo_mon=NaN*zeros(1,240);
  etaglo_aviso_mon=NaN*zeros(1,240);
  etadate=datenum([1992 1 1 12 0 0])+[1:nt]-1;
  for mm=1:240;
    t0=datenum([1992 1+(mm-1) 1 0 0 0]);
    t1=datenum([1992 1+mm 1 0 0 0]);
    tt=find(etadate>t0&etadate<t1);
    etaglo_mon(mm)=mean(etaglo(tt));
    etaglo_aviso_mon(mm)=mean(etaglo_aviso(tt));
  end;
  %
  if doSave; eval(['save ' dirMat 'cost/cost_altimeter_etaglo.mat etaglo*;']); end;
end;

%%%%%%%%%%%%%%%%%
%uncertainties:
%%%%%%%%%%%%%%%%%

sig_mdt=NaN*mygrid.RAC;
file0=[dirSigma nameSigma{1}];
if ~isempty(dir(file0)); sig_mdt=read_bin(file0,1,0); end;
%missing field:
%if useNctiles; sig_mdt=read_nctiles([dirNctiles 'sealevel'],'mdt_sigma'); end;

sig_sladiff_smooth=NaN*mygrid.RAC;
file0=[dirSigma nameSigma{2}];
if ~isempty(dir(file0)); sig_sladiff_smooth=read_bin(file0,1,0); end;
if useNctiles; sig_sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_sigma_r5'); end;

sig_sladiff_point=NaN*mygrid.RAC;
file0=[dirSigma nameSigma{2}];
if ~isempty(dir(file0)); sig_sladiff_point=read_bin(file0,1,0); end;
if useNctiles; sig_sladiff_point=read_nctiles([dirNctiles 'sealevel'],'point_sigma_r5'); end;

sig_mdt(find(sig_mdt==0))=NaN;
sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;

myflds.sig_mdt=sig_mdt;
myflds.sig_sladiff_point=sig_sladiff_point;
myflds.sig_sladiff_smooth=sig_sladiff_smooth;

fprintf('done with init\n');

%%%%%%%%%%%%%%%%%
%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;
        
        %mdt cost function term (misfit plot)
        dif_mdt=NaN*mygrid.RAC;
        file0=[dirEcco 'mdtdiff_smooth.data'];
        if ~isempty(dir(file0)); dif_mdt=read_bin(file0,1,0); end;
        if useNctiles; dif_mdt=read_nctiles([dirNctiles 'sealevel'],'mdt_misfit'); end;
        %
        myflds.dif_mdt=dif_mdt;

        fprintf('done with mdt\n');

     if ~useNctiles;        
        %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;

        tic;
        %pre-load lsc cost function term to add it back to pointwise/1day terms:
        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;
        toc;
     end;%if ~useNctiles;

     if useNctiles;
        if doDifObsOrMod==1; 
            sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_misfit');
        elseif doDifObsOrMod==2;
            sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_sla');
        else;
            sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_misfit');
            sladiff_smooth=sladiff_smooth+read_nctiles([dirNctiles 'sealevel'],'lsc_sla');
        end;
        nRecs=size(sladiff_smooth{1},3);
        TT=1992+[0:nRecs-1]/365.25;
     end;

        fprintf('done with lsc\n');

        %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 ~useNctiles;
            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 back the smoothed values that has been subtracted in cost_gencost_sshv4
            sladiff_smooth_tmp=sladiff_smooth;
            sladiff_smooth_tmp(sladiff_tmp==0)=0;
            sladiff_tmp=sladiff_tmp+sladiff_smooth_tmp;
     end;%if ~useNctiles;

     if useNctiles;
            if doDifObsOrMod==1;
                sladiff_tmp=read_nctiles([dirNctiles 'sealevel'],[myset '_misfit']);
            elseif doDifObsOrMod==2;
                sladiff_tmp=read_nctiles([dirNctiles 'sealevel'],[myset '_sla']);
            else;
                sladiff_tmp=read_nctiles([dirNctiles 'sealevel'],[myset '_misfit']);
                sladiff_tmp=sladiff_tmp+read_nctiles([dirNctiles 'sealevel'],[myset '_sla']);
            end;
            sladiff_tmp(isnan(sladiff_tmp))=0;
     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;

        fprintf('done with point\n');

        %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');

     if ~useNctiles;        
        %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;
     end;%if ~useNctiles;

        %plotting mask : regions of less than 100 observations are omitted
        msk_point(msk_point==0)=NaN;
        msk_100pts=1*(nansum(msk_point,3)>=100);
        msk_100pts(msk_100pts==0)=NaN;
        %store
        myflds.msk_100pts=msk_100pts;

        fprintf('done with msk100\n');
        
     if ~useNctiles;
        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);
        toc;
     end;%if ~useNctiles;

     if useNctiles;
        rms_sladiff_smooth=sqrt(nanmean(sladiff_smooth.^2,3));
        std_sladiff_smooth=nanstd(sladiff_smooth,0,3);
     end;

        %store:
        myflds.rms_sladiff_smooth=rms_sladiff_smooth;
        myflds.std_sladiff_smooth=std_sladiff_smooth;

        fprintf('done with rms\n');
        
     if ~useNctiles;
        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;
        toc;
     end;%if ~useNctiles;

     if useNctiles;
            %store:
            myflds.rms_sladiff_point35d=NaN*rms_sladiff_smooth;
            myflds.std_sladiff_point35d=NaN*std_sladiff_smooth;
     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;

        fprintf(['done with ' suf '\n']);
        
        %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	gforget	1.19	%search for nctiles files
32			useNctiles=0;
33			dirNctiles=[dirModel 'nctiles_remotesensing/sealevel/'];
34			if ~isempty(dir(dirNctiles)); useNctiles=1; end;
35
36			nameSigma={'sigma_MDT_glob_eccollc.bin','slaerr_largescale_r5.err','slaerr_gridscale_r5.err'}; maxLatObs=90
37
38			dirSigma=dirModel;
39			if ~isempty(dir([dirModel 'inputfiles/' nameSigma{1}]));
40	gforget	1.17	dirSigma=[dirModel 'inputfiles/'];
41	gforget	1.14	end;
42
43	gforget	1.19	dirEcco=dirModel;
44			if ~isempty(dir([dirModel 'barfiles']));
45	gforget	1.14	dirEcco=[dirModel 'barfiles' filesep];
46			end;
47
48			if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
49			if isempty(dir(dirMat)); mkdir([dirMat]); end;
50			if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
51			runName=pwd; tmp1=strfind(runName,filesep); runName=runName(tmp1(end)+1:end);
52
53	gforget	1.15	%%%%%%%%%%%%%%%%%
54			%global means:
55			%%%%%%%%%%%%%%%%%
56
57			tmpName=dir([dirEcco 'm_eta_day*data']);
58			avisoName=which('MSL_aviso.bin');
59			if length(tmpName)==1&~isempty(avisoName);
60			m_eta_day=rdmds2gcmfaces([dirEcco tmpName.name(1:end-5)]);
61			avisoName=which('MSL_aviso.bin');
62			etaglo_aviso=read2memory(avisoName)';
63			%global mean
64			nt=size(m_eta_day{1},3); etaglo=NaN*zeros(1,nt);
65			area=mygrid.mskC(:,:,1).*mygrid.RAC; areatot=nansum(area);
66			for tt=1:nt;
67			etaglo(tt)=nansum(area.*m_eta_day(:,:,tt))/areatot;
68			end;
69			%monthly mean
70			etaglo_mon=NaN*zeros(1,240);
71			etaglo_aviso_mon=NaN*zeros(1,240);
72			etadate=datenum([1992 1 1 12 0 0])+[1:nt]-1;
73			for mm=1:240;
74			t0=datenum([1992 1+(mm-1) 1 0 0 0]);
75			t1=datenum([1992 1+mm 1 0 0 0]);
76			tt=find(etadate>t0&etadate<t1);
77			etaglo_mon(mm)=mean(etaglo(tt));
78			etaglo_aviso_mon(mm)=mean(etaglo_aviso(tt));
79			end;
80			%
81			if doSave; eval(['save ' dirMat 'cost/cost_altimeter_etaglo.mat etaglo*;']); end;
82			end;
83
84	gforget	1.14	%%%%%%%%%%%%%%%%%
85	gforget	1.19	%uncertainties:
86			%%%%%%%%%%%%%%%%%
87
88			sig_mdt=NaN*mygrid.RAC;
89			file0=[dirSigma nameSigma{1}];
90			if ~isempty(dir(file0)); sig_mdt=read_bin(file0,1,0); end;
91			%missing field:
92			%if useNctiles; sig_mdt=read_nctiles([dirNctiles 'sealevel'],'mdt_sigma'); end;
93
94			sig_sladiff_smooth=NaN*mygrid.RAC;
95			file0=[dirSigma nameSigma{2}];
96			if ~isempty(dir(file0)); sig_sladiff_smooth=read_bin(file0,1,0); end;
97			if useNctiles; sig_sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_sigma_r5'); end;
98
99			sig_sladiff_point=NaN*mygrid.RAC;
100			file0=[dirSigma nameSigma{2}];
101			if ~isempty(dir(file0)); sig_sladiff_point=read_bin(file0,1,0); end;
102			if useNctiles; sig_sladiff_point=read_nctiles([dirNctiles 'sealevel'],'point_sigma_r5'); end;
103
104			sig_mdt(find(sig_mdt==0))=NaN;
105			sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
106			sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;
107
108			myflds.sig_mdt=sig_mdt;
109			myflds.sig_sladiff_point=sig_sladiff_point;
110			myflds.sig_sladiff_smooth=sig_sladiff_smooth;
111
112	gforget	1.20	fprintf('done with init\n');
113	gforget	1.19
114			%%%%%%%%%%%%%%%%%
115	gforget	1.14	%do computations:
116			%%%%%%%%%%%%%%%%%
117
118	gforget	1.4	for doDifObsOrMod=1:3;
119	gforget	1.14
120	gforget	1.11	if doDifObsOrMod==1; suf='modMobs'; elseif doDifObsOrMod==2; suf='obs'; else; suf='mod'; end;
121
122			if doComp==0;
123			eval(['load ' dirMat 'cost/cost_altimeter_' suf '.mat myflds;']);
124			else;
125
126			%mdt cost function term (misfit plot)
127	gforget	1.19	dif_mdt=NaN*mygrid.RAC;
128			file0=[dirEcco 'mdtdiff_smooth.data'];
129			if ~isempty(dir(file0)); dif_mdt=read_bin(file0,1,0); end;
130			if useNctiles; dif_mdt=read_nctiles([dirNctiles 'sealevel'],'mdt_misfit'); end;
131			%
132	gforget	1.11	myflds.dif_mdt=dif_mdt;
133	gforget	1.19
134	gforget	1.20	fprintf('done with mdt\n');
135	gforget	1.19
136			if ~useNctiles;
137	gforget	1.11	%skip blanks:
138	gforget	1.14	tmp1=dir([dirEcco 'sladiff_smooth*data']);
139	gforget	1.11	nrec=tmp1.bytes/90/1170/4;
140			ttShift=17;
141			%skip 1992
142			listRecs=[1+365-ttShift:17+365*18-ttShift];
143			nRecs=length(listRecs); TT=1993+[0:nRecs-1]/365.25;
144	gforget	1.18
145	gforget	1.19	tic;
146	gforget	1.18	%pre-load lsc cost function term to add it back to pointwise/1day terms:
147			if doDifObsOrMod==1;
148			sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
149			elseif doDifObsOrMod==2;
150			sladiff_smooth=cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
151			else;
152			sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
153			sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
154			end;
155	gforget	1.19	toc;
156			end;%if ~useNctiles;
157
158			if useNctiles;
159			if doDifObsOrMod==1;
160			sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_misfit');
161			elseif doDifObsOrMod==2;
162			sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_sla');
163			else;
164			sladiff_smooth=read_nctiles([dirNctiles 'sealevel'],'lsc_misfit');
165			sladiff_smooth=sladiff_smooth+read_nctiles([dirNctiles 'sealevel'],'lsc_sla');
166			end;
167			nRecs=size(sladiff_smooth{1},3);
168			TT=1992+[0:nRecs-1]/365.25;
169			end;
170
171	gforget	1.20	fprintf('done with lsc\n');
172	gforget	1.19
173	gforget	1.11	%pointwise/1day terms:
174			sladiff_point=repmat(0*mygrid.RAC,[1 1 nRecs]); count_point=sladiff_point;
175			for ii=1:3;
176			if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
177			%topex pointwise misfits:
178	gforget	1.19	if ~useNctiles;
179	gforget	1.11	if doDifObsOrMod==1;
180	gforget	1.14	sladiff_tmp=cost_altimeter_read([dirEcco 'sladiff_' myset '_raw'],ttShift+listRecs);
181	gforget	1.11	elseif doDifObsOrMod==2;
182	gforget	1.14	sladiff_tmp=cost_altimeter_read([dirEcco 'slaobs_' myset '_raw'],ttShift+listRecs);
183	gforget	1.11	else;
184	gforget	1.14	sladiff_tmp=cost_altimeter_read([dirEcco 'sladiff_' myset '_raw'],ttShift+listRecs);
185			sladiff_tmp=sladiff_tmp+cost_altimeter_read([dirEcco 'slaobs_' myset '_raw'],ttShift+listRecs);
186	gforget	1.11	end;
187	gforget	1.19	%add back the smoothed values that has been subtracted in cost_gencost_sshv4
188	gforget	1.18	sladiff_smooth_tmp=sladiff_smooth;
189			sladiff_smooth_tmp(sladiff_tmp==0)=0;
190			sladiff_tmp=sladiff_tmp+sladiff_smooth_tmp;
191	gforget	1.19	end;%if ~useNctiles;
192
193			if useNctiles;
194			if doDifObsOrMod==1;
195			sladiff_tmp=read_nctiles([dirNctiles 'sealevel'],[myset '_misfit']);
196			elseif doDifObsOrMod==2;
197			sladiff_tmp=read_nctiles([dirNctiles 'sealevel'],[myset '_sla']);
198			else;
199			sladiff_tmp=read_nctiles([dirNctiles 'sealevel'],[myset '_misfit']);
200			sladiff_tmp=sladiff_tmp+read_nctiles([dirNctiles 'sealevel'],[myset '_sla']);
201			end;
202			sladiff_tmp(isnan(sladiff_tmp))=0;
203			end;
204
205	gforget	1.11	%add to overall data set:
206			sladiff_point=sladiff_point+sladiff_tmp; count_point=count_point+(sladiff_tmp~=0);
207			%finalize data set:
208			sladiff_tmp(sladiff_tmp==0)=NaN;
209			%compute rms:
210			count_tmp=sum(~isnan(sladiff_tmp),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
211			eval(['myflds.rms_' myset '=sqrt(nansum(sladiff_tmp.^2,3)./count_tmp);']);
212			eval(['myflds.std_' myset '=nanstd(sladiff_tmp,0,3).*msk_tmp;']);
213			end;
214	gforget	1.19
215	gforget	1.20	fprintf('done with point\n');
216	gforget	1.19
217	gforget	1.11	%finalize overall data set:
218			count_point(count_point==0)=NaN;
219			sladiff_point=sladiff_point./count_point;
220			%compute overall rms,std:
221			count_tmp=sum(~isnan(sladiff_point),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
222			myflds.rms_sladiff_point=sqrt(nansum(sladiff_point.^2,3)./count_tmp);
223			myflds.std_sladiff_point=nanstd(sladiff_point,0,3).*msk_tmp;
224			%
225			msk_point=1*(count_point>0);
226			%fill blanks:
227			warning('off','MATLAB:divideByZero');
228			msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>maxLatObs))=NaN;
229			myflds.xtrp_rms_sladiff_point=diffsmooth2D_extrap_inv(myflds.rms_sladiff_point,msk);
230			myflds.xtrp_std_sladiff_point=diffsmooth2D_extrap_inv(myflds.std_sladiff_point,msk);
231			warning('on','MATLAB:divideByZero');
232	gforget	1.19
233			if ~useNctiles;
234	gforget	1.11	%computational mask : only points that were actually observed, and in 35d average
235	gforget	1.14	msk_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
236	gforget	1.11	msk_point35d=1*(msk_point35d~=0);
237			tmp1=sum(msk_point35d==0&msk_point~=0);
238			tmp2=sum(msk_point35d~=0&msk_point~=0);
239			fprintf('after masking : %d omitted, %d retained \n',tmp1,tmp2);
240			msk_point=msk_point.*msk_point35d;
241	gforget	1.19	end;%if ~useNctiles;
242
243			%plotting mask : regions of less than 100 observations are omitted
244	gforget	1.11	msk_point(msk_point==0)=NaN;
245			msk_100pts=1*(nansum(msk_point,3)>=100);
246			msk_100pts(msk_100pts==0)=NaN;
247			%store
248			myflds.msk_100pts=msk_100pts;
249	gforget	1.19
250	gforget	1.20	fprintf('done with msk100\n');
251	gforget	1.11
252	gforget	1.19	if ~useNctiles;
253			tic;
254	gforget	1.11	%lsc cost function term:
255			if doDifObsOrMod==1;
256	gforget	1.14	sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
257	gforget	1.11	elseif doDifObsOrMod==2;
258	gforget	1.14	sladiff_smooth=cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
259	gforget	1.11	else;
260	gforget	1.14	sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
261			sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
262	gforget	1.11	end;
263			%mask missing points:
264			sladiff_smooth(sladiff_smooth==0)=NaN;
265			sladiff_smooth=sladiff_smooth.*msk_point;
266			count_tmp=sum(~isnan(sladiff_smooth),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
267			%compute rms:
268			rms_sladiff_smooth=msk_tmp.*sqrt(nanmean(sladiff_smooth.^2,3));
269			std_sladiff_smooth=msk_tmp.*nanstd(sladiff_smooth,0,3);
270	gforget	1.19	toc;
271			end;%if ~useNctiles;
272
273			if useNctiles;
274			rms_sladiff_smooth=sqrt(nanmean(sladiff_smooth.^2,3));
275			std_sladiff_smooth=nanstd(sladiff_smooth,0,3);
276			end;
277
278	gforget	1.11	%store:
279			myflds.rms_sladiff_smooth=rms_sladiff_smooth;
280			myflds.std_sladiff_smooth=std_sladiff_smooth;
281	gforget	1.19
282	gforget	1.20	fprintf('done with rms\n');
283	gforget	1.11
284	gforget	1.19	if ~useNctiles;
285			tic;
286	gforget	1.11	%pointwise/point35days cost function term:
287			if doDifObsOrMod==1;
288	gforget	1.14	sladiff_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
289	gforget	1.11	elseif doDifObsOrMod==2;
290	gforget	1.14	sladiff_point35d=cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
291	gforget	1.11	else;
292	gforget	1.14	sladiff_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
293			sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
294	gforget	1.11	end;
295			%mask missing points:
296			sladiff_point35d(sladiff_point35d==0)=NaN;
297			sladiff_point35d=sladiff_point35d.*msk_point;
298			count_tmp=sum(~isnan(sladiff_point35d),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
299			%compute rms:
300			rms_sladiff_point35d=msk_tmp.*sqrt(nanmean(sladiff_point35d.^2,3));
301			std_sladiff_point35d=msk_tmp.*nanstd(sladiff_point35d,0,3);
302			%store:
303			myflds.rms_sladiff_point35d=rms_sladiff_point35d;
304			myflds.std_sladiff_point35d=std_sladiff_point35d;
305
306			if 0;
307			%difference between scales
308			rms_sladiff_35dMsmooth=msk_tmp.*sqrt(nanmean((sladiff_point35d-sladiff_smooth).^2,3));
309			std_sladiff_35dMsmooth=msk_tmp.*nanstd((sladiff_point35d-sladiff_smooth),0,3);
310			%store:
311			myflds.rms_sladiff_35dMsmooth=rms_sladiff_35dMsmooth;
312			myflds.std_sladiff_35dMsmooth=std_sladiff_35dMsmooth;
313
314			%difference between scales
315			rms_sladiff_pointMpoint35d=msk_tmp.*sqrt(nanmean((sladiff_point-sladiff_point35d).^2,3));
316			std_sladiff_pointMpoint35d=msk_tmp.*nanstd((sladiff_point-sladiff_point35d),0,3);
317			%store:
318			myflds.rms_sladiff_pointMpoint35d=rms_sladiff_pointMpoint35d;
319			myflds.std_sladiff_pointMpoint35d=std_sladiff_pointMpoint35d;
320			end;
321	gforget	1.19	toc;
322			end;%if ~useNctiles;
323
324			if useNctiles;
325			%store:
326			myflds.rms_sladiff_point35d=NaN*rms_sladiff_smooth;
327			myflds.std_sladiff_point35d=NaN*std_sladiff_smooth;
328			end;
329	gforget	1.11
330			%compute zonal mean/median:
331			for ii=1:4;
332			switch ii;
333			case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
334			case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
335			case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
336			case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
337			end;
338			cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
339			cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
340			end;
341	gforget	1.19
342	gforget	1.20	fprintf(['done with ' suf '\n']);
343	gforget	1.11
344	gforget	1.19	%write to disk:
345	gforget	1.11	if doSave; eval(['save ' dirMat 'cost/cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;
346
347			end;%if doComp
348
349			if doPlot;
350			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];
351			figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
352			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];
353			figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
354			figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
355			figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow;
356			end;
357
358			if doPlot&doPrint;
359			dirFig='../figs/altimeter/'; ff0=gcf-4;
360			for ff=1:4;
361			figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig');
362			eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);
363			eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']);
364			end;
365			end;
366
367	gforget	1.4	end;%for doDifObsOrMod=1:3;
368
369
370
371			function [fldOut]=cost_altimeter_read(fileIn,recIn);
372
373			nrec=length(recIn);
374			global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
375			lrec=siz(1)siz(2)4;
376			myprec='float32';
377
378			fldOut=zeros(siz);
379			fid=fopen([fileIn '.data'],'r','b');
380			for irec=1:nrec;
381	gforget	1.11	status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
382			fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
383	gforget	1.4	end;
384
385			fldOut=convert2gcmfaces(fldOut);
386
387
388