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

Annotation of /MITgcm_contrib/gael/matlab_class/ecco_v4/cost_altimeter_disp.m

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


Revision 1.14 - (hide annotations) (download)
Wed May 11 03:37:11 2016 UTC (9 years, 2 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65w, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
Changes since 1.13: +5 -1 lines
- skip empty plots.

1 gforget 1.5 function []=cost_altimeter_disp(dirMat,choicePlot,suf,dirTex,nameTex);
2     %object: plot the various sea level statistics
3 gforget 1.1 % (std model-obs, model, obs, leading to cost function terms)
4 gforget 1.4 %inputs: dirMat is the model run directory
5 gforget 1.5 % choicePlot is 1 (rms) 2 (prior uncertainty) or 3 (cost)
6     % suf is 'modMobs', 'obs' or 'mod'
7     %optional: dirTex is the directory where tex and figures files are created
8     % (if not specified then display all results to screen instead)
9     % nameTex is the tex file name (default : 'myPlots')
10 gforget 1.1
11 gforget 1.5 gcmfaces_global;
12    
13     %backward compatibility test
14     test1=~isempty(dir([dirMat 'basic_diags_ecco_mygrid.mat']));
15     test2=~isempty(dir([dirMat 'diags_grid_parms.mat']));
16     if ~test1&~test2;
17 gforget 1.8 error('missing diags_grid_parms.mat')
18 gforget 1.5 elseif test2;
19 gforget 1.8 nameGrid='diags_grid_parms.mat';
20 gforget 1.5 else;
21 gforget 1.8 nameGrid='basic_diags_ecco_mygrid.mat';
22 gforget 1.5 end;
23 gforget 1.1
24 gforget 1.4 %here we always reload the grid from dirMat to make sure the same one is used throughout
25 gforget 1.5 eval(['load ' dirMat nameGrid ';']);
26 gforget 1.1
27 gforget 1.4 %determine if and where to create tex and figures files
28     dirMat=[dirMat '/'];
29     if isempty(who('dirTex'));
30 gforget 1.8 addToTex=0;
31 gforget 1.4 else;
32 gforget 1.8 if ~ischar(dirTex); error('mis-specified dirTex'); end;
33     addToTex=1;
34     if isempty(who('nameTex')); nameTex='myPlots'; end;
35 gforget 1.11 fileTex=[dirTex filesep nameTex '.tex'];
36 gforget 1.4 end;
37 gforget 1.1
38     %%%%%%%%%%%%%%%
39     %define pathes:
40     %%%%%%%%%%%%%%%
41    
42 gforget 1.4 if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
43 gforget 1.10 runName=pwd; tmp1=strfind(runName,filesep); runName=runName(tmp1(end)+1:end);
44 gforget 1.1
45     %%%%%%%%%%%%%%%%%
46     %do computations:
47     %%%%%%%%%%%%%%%%%
48    
49 gforget 1.8 if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
50    
51 gforget 1.12 %global means
52     if choicePlot==0;
53     nmMat='diags_set_C';
54     test0=~isempty(dir([dirMat 'cost_altimeter_etaglo.mat']));
55     test0=test0&~isempty(dir([dirMat '../' nmMat filesep nmMat '_*.mat']));
56     if test0;
57     eval(['load ' dirMat 'cost_altimeter_etaglo.mat;']);
58     etaglo_offset=nanmedian(etaglo_aviso_mon-etaglo_mon);
59     alldiag=diags_read_from_mat([dirMat '../' nmMat filesep],[nmMat '_*.mat']);
60     alldiag=diags_read_from_mat([dirMat '../' nmMat filesep],[nmMat '_*.mat'],'SSHglo');
61     %
62     tt=alldiag.listTimes;
63     etaglo_mass_mon=alldiag.SSHglo;
64 gforget 1.13 %
65     if 0;
66     %
67     ii=max(find(tt<2005));
68     etaglo_mon=etaglo_mon(ii+1:end);
69     etaglo_mass_mon=etaglo_mass_mon(ii+1:end);
70     etaglo_aviso_mon=etaglo_aviso_mon(ii+1:end);
71     etaglo_offset=0; tt=tt(ii+1:end);
72     %
73     etaglo_cy=zeros(1,12); etaglo_mass_cy=zeros(1,12); etaglo_aviso_cy=zeros(1,12);
74     for jj=1:12;
75     etaglo_cy(jj)=mean(etaglo_mon(jj:12:end));
76     etaglo_aviso_cy(jj)=mean(etaglo_aviso_mon(jj:12:end));
77     etaglo_mass_cy(jj)=mean(etaglo_mass_mon(jj:12:end));
78     end;
79     etaglo_cy=repmat(etaglo_cy,[1 8]); etaglo_mon=etaglo_mon-etaglo_cy(1:85);
80     etaglo_aviso_cy=repmat(etaglo_aviso_cy,[1 8]); etaglo_aviso_mon=etaglo_aviso_mon-etaglo_aviso_cy(1:85);
81     etaglo_mass_cy=repmat(etaglo_mass_cy,[1 8]); etaglo_mass_mon=etaglo_mass_mon-etaglo_mass_cy(1:85);
82     %
83     etaglo_mon=etaglo_mon-etaglo_mon(1);
84     etaglo_mass_mon=etaglo_mass_mon-etaglo_mass_mon(1);
85     etaglo_aviso_mon=etaglo_aviso_mon-etaglo_aviso_mon(1);
86     %
87     end;%if 0;
88     %
89 gforget 1.12 figure; plot(tt,etaglo_mon+etaglo_offset,'k','LineWidth',1); hold on;
90     plot(tt,etaglo_mass_mon+etaglo_offset,'b','LineWidth',1);
91     plot(tt,etaglo_mon-etaglo_mass_mon,'r','LineWidth',1);
92     plot(tt,etaglo_aviso_mon,'m','LineWidth',1);
93     legend('sea level','mass term','steric term','aviso','Location','NorthWest');
94     myCaption={'global sea level (m)'};
95     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
96     end;
97     %skip the rest
98     return;
99     end;
100    
101     %load results:
102 gforget 1.4 eval(['load ' dirMat 'cost_altimeter_' suf '.mat myflds;']);
103 gforget 1.1
104 gforget 1.8 %mask for plotting
105     if isfield(myflds,'msk_100pts');
106     MSK=myflds.msk_100pts;
107     else;
108     MSK=mygrid.mskC(:,:,1);
109     end;
110     %hack to omit any masking : MSK(:)=1;
111    
112 gforget 1.1 if strcmp(suf,'modMobs'); tit='modeled-observed';
113     elseif strcmp(suf,'obs'); tit='observed';
114     elseif strcmp(suf,'mod'); tit='modeled';
115     else; error('unknown field');
116     end
117    
118 gforget 1.9 if choicePlot==1; tit=[tit ' log(variance)']; uni='(m$^2$)';
119     elseif choicePlot==2; tit=' log(prior error variance)'; uni='(m$^2$)';
120 gforget 1.8 else; tit=[tit ' cost']; uni='';
121 gforget 1.1 end;
122    
123     if choicePlot==1;%rms
124 gforget 1.8
125     if strcmp(suf,'modMobs');
126     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];
127     figure; m_map_gcmfaces(100*MSK.*myflds.dif_mdt,0,{'myCaxis',100*cc}); drawnow;
128 gforget 1.9 myCaption={'mean dynamic topography misfit (cm)'};
129 gforget 1.8 if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
130     end;
131    
132     cc=[-4:0.2:-1]-0.5;
133     figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_smooth.^2)),0,{'myCaxis',cc}); drawnow;
134     myCaption={tit,'-- sea level anomaly ',uni,' -- large space/time scales'};
135     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
136    
137     cc=[-4:0.2:-1];
138     figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_point.^2)),0,{'myCaxis',cc}); drawnow;
139     myCaption={tit,'-- sea level anomaly ',uni,' -- pointwise'};
140     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
141    
142     if 0;
143     figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_point35d.^2)),0,{'myCaxis',cc}); drawnow;
144     myCaption={tit,'-- sea level anomaly ',uni,' -- large time scales'};
145     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
146    
147     figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_35dMsmooth.^2)),0,{'myCaxis',cc}); drawnow;
148     myCaption={tit,'-- sea level anomaly ',uni,' -- point35d minus lsc'};
149     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
150    
151     figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_pointMpoint35d.^2)),0,{'myCaxis',cc}); drawnow;
152     myCaption={tit,'-- sea level anomaly ',uni,' -- point minus point35d'};
153     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
154     end;
155    
156 gforget 1.1 elseif choicePlot==2;%uncertainty fields
157 gforget 1.14
158     if sum(~isnan(myflds.sig_mdt))>0;
159 gforget 1.8 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];
160     figure; m_map_gcmfaces(100*myflds.sig_mdt,0,{'myCaxis',100*cc}); drawnow;
161 gforget 1.9 myCaption={'mean dynamic topography prior uncertainty (cm)'};
162 gforget 1.8 if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
163 gforget 1.14 end;
164 gforget 1.8
165     cc=[-4:0.2:-1]-0.5;
166     figure; m_map_gcmfaces(log10(MSK.*(myflds.sig_sladiff_smooth.^2)),0,{'myCaxis',cc}); drawnow;
167     myCaption={tit,'-- sea level anomaly ',uni,' -- large space/time scales'};
168     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
169    
170     cc=[-4:0.2:-1];
171     figure; m_map_gcmfaces(log10(MSK.*(myflds.sig_sladiff_point.^2)),0,{'myCaxis',cc}); drawnow;
172     myCaption={tit,'-- sea level anomaly ',uni,' -- pointwise'};
173     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
174    
175 gforget 1.1 else;%cost
176 gforget 1.8
177     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]*100;
178    
179 gforget 1.14 if sum(~isnan(myflds.sig_mdt))>0;
180 gforget 1.8 figure; m_map_gcmfaces(MSK.*((myflds.dif_mdt.^2)./(myflds.sig_mdt.^2)),0,{'myCaxis',cc}); drawnow;
181     myCaption={tit,'-- mean dynamic topography ',uni};
182     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
183 gforget 1.14 end;
184 gforget 1.8
185     figure; m_map_gcmfaces(MSK.*((myflds.rms_sladiff_smooth.^2)./(myflds.sig_sladiff_smooth.^2)),0,{'myCaxis',cc}); drawnow;
186     myCaption={tit,'-- sea level anomaly ',uni,' -- large space/time scales'};
187     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
188    
189     figure; m_map_gcmfaces(MSK.*((myflds.rms_sladiff_point.^2)./(myflds.sig_sladiff_point.^2)),0,{'myCaxis',cc}); drawnow;
190     myCaption={tit,'-- sea level anomaly ',uni,' -- pointwise'};
191     if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
192    
193 gforget 1.1 end;
194    
195    
196    
197    
198    

  ViewVC Help
Powered by ViewVC 1.1.22