/[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.4 - (hide annotations) (download)
Fri Jul 1 16:07:26 2011 UTC (14 years ago) by gforget
Branch: MAIN
Changes since 1.3: +93 -33 lines 
ecco_v4/basic_diags_ecco.m
        make a function of dirModel,lChunk,iChunkm[,reLoadGrid]
        break records loop in chunks, and read nChnunk records at once
        write output to dirModel/mat/sic_diags_ecco_*.mat
        get grid from GRID/ if needed or as specified by reLoadGrid
ecco_v4/basic_diags_ecco_disp.m
        now a function of dirModel,addToTex
        dont reload grid
        m_map_gcmfaces updates
        added write2tex calls
ecco_v4/cost_altimeter.m
        function of dirModel
        added doDifObsOrMod loop
        get grid from GRID/ if needed
        the display part is kept for now, but the it will be
                replaced with new cost_altimeter_disp.m
        include faster cost_altimeter_read
        write output to dirModel/mat/cost_altimeter_*.mat
added functions
        ecco_v4/cost_altimeter_disp.m
        ecco_v4/cost_bp.m
        ecco_v4/cost_sst.m
        ecco_v4/plot_driver.m
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     %store:
51     myflds.dif_mdt=dif_mdt;
52     myflds.sig_mdt=sig_mdt;
53    
54 gforget 1.4 %skip blanks:
55     listRecs=[1:7:6147];
56    
57 gforget 1.1 toc; tic;
58     %lsc cost function term:
59 gforget 1.4 if doDifObsOrMod==1;
60     sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
61     elseif doDifObsOrMod==2;
62     sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
63     else;
64     sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
65     sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
66     end;
67     %mask missing points:
68     sladiff_smooth(sladiff_smooth==0)=NaN;
69 gforget 1.1 %compute rms:
70     rms_sladiff_smooth=sqrt(nanmean(sladiff_smooth.^2,3));
71     %get weight:
72     sig_sladiff_smooth=v4_read_bin([dirData 'sigma_SLA_smooth_eccollc.bin'],1,0);
73     %store:
74     myflds.rms_sladiff_smooth=rms_sladiff_smooth;
75     myflds.sig_sladiff_smooth=sig_sladiff_smooth;
76    
77     toc; tic;
78 gforget 1.4 %pointwise/point35days cost function term:
79     if doDifObsOrMod==1;
80     sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
81     elseif doDifObsOrMod==2;
82     sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
83     else;
84     sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
85     sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
86     end;
87     %mask missing points:
88     sladiff_point35d(sladiff_point35d==0)=NaN;
89 gforget 1.1 %compute rms:
90 gforget 1.4 rms_sladiff_point35d=sqrt(nanmean(sladiff_point35d.^2,3));
91 gforget 1.1 %store:
92 gforget 1.4 myflds.rms_sladiff_point35d=rms_sladiff_point35d;
93 gforget 1.1
94     toc; tic;
95     %pointwise/1day terms:
96     sum_all=0*mygrid.XC; msk_all=sum_all;
97     for ii=1:3;
98     if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
99     %topex pointwise misfits:
100 gforget 1.4 if doDifObsOrMod==1;
101     sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
102     elseif doDifObsOrMod==2;
103     sladiff_point=cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
104     else;
105     sladiff_point=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],listRecs);
106     sladiff_point=sladiff_point+cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],listRecs);
107     end;
108 gforget 1.1 %compute rms:
109     msk_tmp=1*(sladiff_point~=0);
110     msk_tmp=sum(msk_tmp,3); sum_tmp=sum(sladiff_point.^2,3);
111     sum_all=sum_all+sum_tmp; msk_all=msk_all+msk_tmp;
112     eval(['rms_' myset '=sum_tmp./msk_tmp;']);
113     end;
114     %compute overall rms:
115     rms_sladiff_point=sqrt(sum_all./msk_all);
116     %fill blanks:
117     rms_sladiff_point=diffsmooth2D_extrap_inv(rms_sladiff_point,mygrid.mskC(:,:,1));
118     %get weight:
119     sig_sladiff_point=v4_read_bin([dirData 'sigma_SLA_PWS07r2_glob_eccollc.bin'],1,0);
120     %store:
121     myflds.rms_tp=rms_tp; myflds.rms_gfo=rms_gfo; myflds.rms_ers=rms_ers;
122     myflds.rms_sladiff_point=rms_sladiff_point;
123     myflds.sig_sladiff_point=sig_sladiff_point;
124    
125     %compute zonal mean/median:
126 gforget 1.4 for ii=1:4;
127 gforget 1.1 switch ii;
128     case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
129     case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
130 gforget 1.4 case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
131     case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
132 gforget 1.1 end;
133 gforget 1.3 cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
134     cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
135 gforget 1.1 end;
136    
137     toc; %write to disk:
138 gforget 1.4 if doSave; eval(['save ' dirOut 'cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;
139 gforget 1.1
140     end;%if doComp
141    
142     if doPlot;
143 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];
144     figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
145     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];
146     figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
147     figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
148     figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow;
149 gforget 1.1 end;
150    
151     if doPlot&doPrint;
152     dirFig='../figs/altimeter/'; ff0=gcf-4;
153     for ff=1:4;
154 gforget 1.4 figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig');
155     eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);
156     eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']);
157 gforget 1.1 end;
158     end;
159    
160 gforget 1.4 end;%for doDifObsOrMod=1:3;
161    
162    
163    
164     function [fldOut]=cost_altimeter_read(fileIn,recIn);
165    
166     nrec=length(recIn);
167     global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
168     lrec=siz(1)*siz(2)*4;
169     myprec='float32';
170    
171     fldOut=zeros(siz);
172     fid=fopen([fileIn '.data'],'r','b');
173     for irec=1:nrec;
174     status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
175     fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
176     end;
177    
178     fldOut=convert2gcmfaces(fldOut);
179    
180    
181    
  ViewVC Help
Powered by ViewVC 1.1.22