/[MITgcm]/MITgcm_contrib/gael/matlab_class/ecco_v4/cost_altimeter.m
ViewVC logotype

Contents of /MITgcm_contrib/gael/matlab_class/ecco_v4/cost_altimeter.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.21 - (show annotations) (download)
Fri Aug 5 19:29:40 2016 UTC (8 years, 11 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
Changes since 1.20: +1 -1 lines
- bug fix (reported by A. Nguyen).

1 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
9 doComp=1;
10 if doComp==1;
11 doSave=1;
12 doPlot=0;
13 else;
14 doSave=0;
15 doPlot=1;
16 end;
17 doPrint=0;
18
19 %%%%%%%%%%%
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 %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 dirSigma=[dirModel 'inputfiles/'];
41 end;
42
43 dirEcco=dirModel;
44 if ~isempty(dir([dirModel 'barfiles']));
45 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 %%%%%%%%%%%%%%%%%
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 %%%%%%%%%%%%%%%%%
85 %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{3}];
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 fprintf('done with init\n');
113
114 %%%%%%%%%%%%%%%%%
115 %do computations:
116 %%%%%%%%%%%%%%%%%
117
118 for doDifObsOrMod=1:3;
119
120 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 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 myflds.dif_mdt=dif_mdt;
133
134 fprintf('done with mdt\n');
135
136 if ~useNctiles;
137 %skip blanks:
138 tmp1=dir([dirEcco 'sladiff_smooth*data']);
139 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
145 tic;
146 %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 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 fprintf('done with lsc\n');
172
173 %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 if ~useNctiles;
179 if doDifObsOrMod==1;
180 sladiff_tmp=cost_altimeter_read([dirEcco 'sladiff_' myset '_raw'],ttShift+listRecs);
181 elseif doDifObsOrMod==2;
182 sladiff_tmp=cost_altimeter_read([dirEcco 'slaobs_' myset '_raw'],ttShift+listRecs);
183 else;
184 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 end;
187 %add back the smoothed values that has been subtracted in cost_gencost_sshv4
188 sladiff_smooth_tmp=sladiff_smooth;
189 sladiff_smooth_tmp(sladiff_tmp==0)=0;
190 sladiff_tmp=sladiff_tmp+sladiff_smooth_tmp;
191 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 %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
215 fprintf('done with point\n');
216
217 %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
233 if ~useNctiles;
234 %computational mask : only points that were actually observed, and in 35d average
235 msk_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
236 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 end;%if ~useNctiles;
242
243 %plotting mask : regions of less than 100 observations are omitted
244 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
250 fprintf('done with msk100\n');
251
252 if ~useNctiles;
253 tic;
254 %lsc cost function term:
255 if doDifObsOrMod==1;
256 sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
257 elseif doDifObsOrMod==2;
258 sladiff_smooth=cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
259 else;
260 sladiff_smooth=cost_altimeter_read([dirEcco 'sladiff_smooth'],listRecs);
261 sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirEcco 'slaobs_smooth'],listRecs);
262 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 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 %store:
279 myflds.rms_sladiff_smooth=rms_sladiff_smooth;
280 myflds.std_sladiff_smooth=std_sladiff_smooth;
281
282 fprintf('done with rms\n');
283
284 if ~useNctiles;
285 tic;
286 %pointwise/point35days cost function term:
287 if doDifObsOrMod==1;
288 sladiff_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
289 elseif doDifObsOrMod==2;
290 sladiff_point35d=cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
291 else;
292 sladiff_point35d=cost_altimeter_read([dirEcco 'sladiff_raw'],listRecs);
293 sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirEcco 'slaobs_raw'],listRecs);
294 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 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
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
342 fprintf(['done with ' suf '\n']);
343
344 %write to disk:
345 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 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 status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
382 fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
383 end;
384
385 fldOut=convert2gcmfaces(fldOut);
386
387
388

  ViewVC Help
Powered by ViewVC 1.1.22