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

Annotation 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.5 - (hide annotations) (download)
Fri Jul 8 13:18:04 2011 UTC (14 years ago) by gforget
Branch: MAIN
Changes since 1.4: +10 -2 lines 
- added: cost_xx.m that compute/plot the mean and variance of controls.
- cost_bp.m, cost_sst.m, cost_altimeter_disp.m : revised captions for wrote2tex.
- plot_driver.m : add title, sections, cost_xx. split basic_diags_ecco_disp.
- cost_altimeter.m : mask out data poleward of 66. resolve divide by zero
  warnings by applying masks, and switching off the warning.
- basic_diags_ecco.m : split in three using new setDiags paramater; speed-up I/O
  by using rdmds2workspace_list.m; fix heat and salt transport computations (missing
  factors); include time averaging for setDiags='B'; save mygrid to basic_diags_ecco_mygrid.mat
- basic_diags_ecco_disp.m : split in three using new setDiags paramater; load mygrid
  from basic_diags_ecco_mygrid.mat; improve captions for write2tex; set time averaging
  period to 92-07; include computation of zonal mean climatology for setDiags='C';
  use depthStretchPlot for section; add air-sea fluxes and ice plots; fix transports
  time series plots;
1 gforget 1.4 function []=cost_altimeter(dirModel);
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    
6     % doComp is the switch from computation to plot
7     doComp=1;
8     if doComp==1;
9     doSave=1;
10     doPlot=0;
11     else;
12     doSave=0;
13     doPlot=1;
14     end;
15     doPrint=0;
16    
17     for doDifObsOrMod=1:3;
18 gforget 1.1
19 gforget 1.4 if doDifObsOrMod==1; suf='modMobs'; elseif doDifObsOrMod==2; suf='obs'; else; suf='mod'; end;
20 gforget 1.1
21     %%%%%%%%%%%
22     %load grid:
23     %%%%%%%%%%%
24    
25 gforget 1.4 dirGrid='GRID/';
26 gforget 1.1 global mygrid;
27 gforget 1.4 gcmfaces_path; grid_load(dirGrid,5,'compact');
28     gcmfaces_lines_zonal;
29 gforget 1.1
30     %%%%%%%%%%%%%%%
31     %define pathes:
32     %%%%%%%%%%%%%%%
33    
34     dirData='/net/altix3700/raid4/gforget/mysetups/ecco_v4/INPUTS_90x50/';
35 gforget 1.4 dirOut=[dirModel 'mat/']; if isempty(dir(dirOut)); eval(['mkdir ' dirOut ';']); end;
36 gforget 1.1 runName=pwd; tmp1=strfind(runName,'/'); runName=runName(tmp1(end)+1:end);
37    
38     %%%%%%%%%%%%%%%%%
39     %do computations:
40     %%%%%%%%%%%%%%%%%
41    
42     if doComp==0;
43 gforget 1.4 eval(['load ' dirOut 'cost_altimeter_' suf '.mat myflds;']);
44 gforget 1.1 else;
45    
46     tic;
47     %mdt cost function term (misfit plot)
48 gforget 1.4 dif_mdt=rdmds2gcmfaces([dirModel 'barfiles/mdtdiff_smooth']);
49 gforget 1.1 sig_mdt=v4_read_bin([dirData 'sigma_MDT_glob_eccollc.bin'],1,0);
50 gforget 1.5 sig_mdt(find(sig_mdt==0))=NaN;
51 gforget 1.1 %store:
52     myflds.dif_mdt=dif_mdt;
53     myflds.sig_mdt=sig_mdt;
54    
55 gforget 1.4 %skip blanks:
56     listRecs=[1:7:6147];
57    
58 gforget 1.1 toc; tic;
59     %lsc cost function term:
60 gforget 1.4 if doDifObsOrMod==1;
61     sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
62     elseif doDifObsOrMod==2;
63     sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
64     else;
65     sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
66     sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
67     end;
68     %mask missing points:
69     sladiff_smooth(sladiff_smooth==0)=NaN;
70 gforget 1.1 %compute rms:
71     rms_sladiff_smooth=sqrt(nanmean(sladiff_smooth.^2,3));
72     %get weight:
73     sig_sladiff_smooth=v4_read_bin([dirData 'sigma_SLA_smooth_eccollc.bin'],1,0);
74 gforget 1.5 sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;
75 gforget 1.1 %store:
76     myflds.rms_sladiff_smooth=rms_sladiff_smooth;
77     myflds.sig_sladiff_smooth=sig_sladiff_smooth;
78    
79     toc; tic;
80 gforget 1.4 %pointwise/point35days cost function term:
81     if doDifObsOrMod==1;
82     sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
83     elseif doDifObsOrMod==2;
84     sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
85     else;
86     sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
87     sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
88     end;
89     %mask missing points:
90     sladiff_point35d(sladiff_point35d==0)=NaN;
91 gforget 1.1 %compute rms:
92 gforget 1.4 rms_sladiff_point35d=sqrt(nanmean(sladiff_point35d.^2,3));
93 gforget 1.1 %store:
94 gforget 1.4 myflds.rms_sladiff_point35d=rms_sladiff_point35d;
95 gforget 1.1
96     toc; tic;
97     %pointwise/1day terms:
98     sum_all=0*mygrid.XC; msk_all=sum_all;
99     for ii=1:3;
100     if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
101     %topex pointwise misfits:
102 gforget 1.4 if doDifObsOrMod==1;
103     sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
104     elseif doDifObsOrMod==2;
105     sladiff_point=cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
106     else;
107     sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
108     sladiff_point=sladiff_point+cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
109     end;
110 gforget 1.1 %compute rms:
111     msk_tmp=1*(sladiff_point~=0);
112     msk_tmp=sum(msk_tmp,3); sum_tmp=sum(sladiff_point.^2,3);
113 gforget 1.5 sum_all=sum_all+sum_tmp; msk_all=msk_all+msk_tmp;
114     msk_tmp(find(msk_tmp==0))=NaN;
115 gforget 1.1 eval(['rms_' myset '=sum_tmp./msk_tmp;']);
116     end;
117     %compute overall rms:
118 gforget 1.5 msk_all(find(msk_all==0))=NaN;
119 gforget 1.1 rms_sladiff_point=sqrt(sum_all./msk_all);
120     %fill blanks:
121 gforget 1.5 warning('off','MATLAB:divideByZero');
122     msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>66))=NaN;
123     rms_sladiff_point=diffsmooth2D_extrap_inv(rms_sladiff_point,msk);
124     warning('on','MATLAB:divideByZero');
125 gforget 1.1 %get weight:
126     sig_sladiff_point=v4_read_bin([dirData 'sigma_SLA_PWS07r2_glob_eccollc.bin'],1,0);
127 gforget 1.5 sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
128 gforget 1.1 %store:
129     myflds.rms_tp=rms_tp; myflds.rms_gfo=rms_gfo; myflds.rms_ers=rms_ers;
130     myflds.rms_sladiff_point=rms_sladiff_point;
131     myflds.sig_sladiff_point=sig_sladiff_point;
132    
133     %compute zonal mean/median:
134 gforget 1.4 for ii=1:4;
135 gforget 1.1 switch ii;
136     case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
137     case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
138 gforget 1.4 case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
139     case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
140 gforget 1.1 end;
141 gforget 1.3 cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
142     cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
143 gforget 1.1 end;
144    
145     toc; %write to disk:
146 gforget 1.4 if doSave; eval(['save ' dirOut 'cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;
147 gforget 1.1
148     end;%if doComp
149    
150     if doPlot;
151 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];
152     figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
153     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];
154     figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
155     figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
156     figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow;
157 gforget 1.1 end;
158    
159     if doPlot&doPrint;
160     dirFig='../figs/altimeter/'; ff0=gcf-4;
161     for ff=1:4;
162 gforget 1.4 figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig');
163     eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);
164     eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']);
165 gforget 1.1 end;
166     end;
167    
168 gforget 1.4 end;%for doDifObsOrMod=1:3;
169    
170    
171    
172     function [fldOut]=cost_altimeter_read(fileIn,recIn);
173    
174     nrec=length(recIn);
175     global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
176     lrec=siz(1)*siz(2)*4;
177     myprec='float32';
178    
179     fldOut=zeros(siz);
180     fid=fopen([fileIn '.data'],'r','b');
181     for irec=1:nrec;
182     status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
183     fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
184     end;
185    
186     fldOut=convert2gcmfaces(fldOut);
187    
188    
189    
  ViewVC Help
Powered by ViewVC 1.1.22