1 | function stats = evaluate_ckd_lw_fluxes(ref_file, ckd_file_in, model_name, scenario_title) |
---|
2 | % Evaluate the fluxes and heating rates from a CKD model with |
---|
3 | % line-by-line benchmark calculations for a CKDMIP scenario |
---|
4 | |
---|
5 | ickd_stats = 1; % Which CKD model do we show the stats for (0=none) |
---|
6 | |
---|
7 | if iscell(ckd_file_in) |
---|
8 | ckd_file = ckd_file_in; |
---|
9 | else |
---|
10 | ckd_file = {ckd_file_in}; |
---|
11 | end |
---|
12 | |
---|
13 | if iscell(model_name) |
---|
14 | titles = model_name; |
---|
15 | else |
---|
16 | titles = {model_name} |
---|
17 | end |
---|
18 | |
---|
19 | combined_titles = ''; |
---|
20 | for ie = 1:length(titles) |
---|
21 | combined_titles = [combined_titles titles{ie} '_']; |
---|
22 | end |
---|
23 | ickd = 1:length(ckd_file); |
---|
24 | |
---|
25 | col_ref = 'k'; |
---|
26 | cols = {'r','b','g','m','c'}; |
---|
27 | alpha= [0.2 0.3 0.25 0.25 0.1 0.1]; |
---|
28 | |
---|
29 | dflux_axis = [-4 4]; |
---|
30 | %dfluxscatter_axis = [-8 8]; |
---|
31 | dfluxscatter_axis = [-4 4]; |
---|
32 | hr_axis = [-1 1]; |
---|
33 | |
---|
34 | if length(ckd_file) > 1 |
---|
35 | do_legend = 1; |
---|
36 | else |
---|
37 | do_legend = 0; |
---|
38 | end |
---|
39 | |
---|
40 | legfontsize = 7; |
---|
41 | error_save = []; |
---|
42 | |
---|
43 | |
---|
44 | [ref,attr_ref] = loadnc(ref_file); |
---|
45 | if ~exist('scenario_title','var') |
---|
46 | scenario_title = attr_ref.global.scenario; |
---|
47 | end |
---|
48 | |
---|
49 | hr_ref = calc_hr(ref,'lw')'; |
---|
50 | |
---|
51 | for ie = 1:length(ckd_file) |
---|
52 | ckd{ie} = loadnc(ckd_file{ie}); |
---|
53 | hr_ckd{ie} = calc_hr(ckd{ie},'lw')'; |
---|
54 | end |
---|
55 | |
---|
56 | % Pressure in hPa at full levels and half levels |
---|
57 | p_fl_hPa = 0.01.*0.5.*(ref.pressure_hl(1:end-1,:) + ref.pressure_hl(2:end,:)); |
---|
58 | p_hl_hPa = 0.01.*ref.pressure_hl; |
---|
59 | |
---|
60 | iprof = 1:size(ref.pressure_hl,2); |
---|
61 | |
---|
62 | xloc = 0.025;yloc = 0.95; |
---|
63 | |
---|
64 | if nargout == 0 |
---|
65 | do_plot = 1; |
---|
66 | else |
---|
67 | do_plot = 0; |
---|
68 | end |
---|
69 | |
---|
70 | if do_plot |
---|
71 | clf |
---|
72 | set(gcf,'paperposition',[0.5 0.5 27 25]); |
---|
73 | set(gcf,'defaultaxeslayer','top','defaultlinelinewidth',0.5,'defaulttextfontsize',10); |
---|
74 | |
---|
75 | % Plot reference flux and heating-rate profiles |
---|
76 | subplot(3,3,1) |
---|
77 | semilogy(ref.flux_up_lw(:,iprof), p_hl_hPa(:,iprof), col_ref); |
---|
78 | set(gca,'ydir','reverse'); |
---|
79 | ylim([0.01 1000]); |
---|
80 | ylabel('Pressure (hPa)'); |
---|
81 | xlabel('Upwelling longwave flux (W m^{-2})'); |
---|
82 | title('Reference profiles') |
---|
83 | text(xloc, yloc, '\bf(a)','units','normalized'); |
---|
84 | |
---|
85 | subplot(3,3,4) |
---|
86 | semilogy(ref.flux_dn_lw(:,iprof), p_hl_hPa, col_ref); |
---|
87 | set(gca,'ydir','reverse'); |
---|
88 | ylim([0.01 1000]); |
---|
89 | ylabel('Pressure (hPa)'); |
---|
90 | xlabel('Downwelling longwave flux (W m^{-2})'); |
---|
91 | text(xloc, yloc, '\bf(d)','units','normalized'); |
---|
92 | |
---|
93 | subplot(3,3,7) |
---|
94 | semilogy(hr_ref(:,iprof), p_fl_hPa(:,iprof), col_ref); |
---|
95 | set(gca,'ydir','reverse'); |
---|
96 | ylim([0.01 1000]); |
---|
97 | ylabel('Pressure (hPa)'); |
---|
98 | xlabel('Heating rate (K d^{-1})'); |
---|
99 | xlim([-25 5]); |
---|
100 | text(xloc, yloc, '\bf(g)','units','normalized'); |
---|
101 | |
---|
102 | % Plot random errors and biases in flux and heating-rate profiles |
---|
103 | subplot(3,3,8) |
---|
104 | end |
---|
105 | |
---|
106 | ileg = 1; |
---|
107 | leg = {}; |
---|
108 | for ie = ickd |
---|
109 | err = calc_hr_error(p_hl_hPa(:,iprof), hr_ckd{ie}(:,iprof),hr_ref(:,iprof), [0.02 4]); |
---|
110 | err4 = calc_hr_error(p_hl_hPa(:,iprof), hr_ckd{ie}(:,iprof),hr_ref(:,iprof), [4 1100]); |
---|
111 | leg{ileg} = [titles{ie} ' (RMSE=' num2str(err,'%0.3f') ' K d^{-1})']; |
---|
112 | ileg = ileg+1; |
---|
113 | error_save(1,ie) = err; |
---|
114 | error_save(2,ie) = err4; |
---|
115 | |
---|
116 | if do_plot |
---|
117 | plot(-1,-1,cols{ie},'linewidth',1.5); |
---|
118 | hold on |
---|
119 | end |
---|
120 | end |
---|
121 | |
---|
122 | |
---|
123 | if do_plot |
---|
124 | for ie = ickd |
---|
125 | hr_error = hr_ckd{ie}(:,iprof) - hr_ref(:,iprof); |
---|
126 | hr_bias = mean(hr_error'); |
---|
127 | hr_ci = std(hr_error').*1.96; |
---|
128 | hr_errmin = min(hr_error,[],2); |
---|
129 | hr_errmax = max(hr_error,[],2); |
---|
130 | pp = mean(p_fl_hPa'); |
---|
131 | h=fill([hr_bias+hr_ci flip(hr_bias-hr_ci)],... |
---|
132 | [pp flip(pp)],'r'); |
---|
133 | set(h,'facecolor',cols{ie}(1),'edgecolor','none','facealpha',alpha(ie)); |
---|
134 | hold on |
---|
135 | plot(hr_bias,pp,cols{ie},'linewidth',1.5); |
---|
136 | % plot(hr_errmin,pp,[cols{ie}(1) '--']) |
---|
137 | % plot(hr_errmax,pp,[cols{ie}(1) '--']) |
---|
138 | end |
---|
139 | plot(hr_axis,[0.02 0.02],'k:'); |
---|
140 | plot(hr_axis,[4 4],'k:'); |
---|
141 | set(gca,'ydir','reverse'); |
---|
142 | set(gca,'yscale','log'); |
---|
143 | ylim([0.01 1000]); |
---|
144 | ylabel('Pressure (hPa)'); |
---|
145 | xlabel('Heating rate error (K d^{-1})'); |
---|
146 | xlim(hr_axis); |
---|
147 | plot([0 0],[0.01 1000],'k:'); |
---|
148 | %if do_legend |
---|
149 | % set(legend(leg,'location','south'),'fontsize',legfontsize','box','on'); |
---|
150 | %end |
---|
151 | text(xloc, yloc, '\bf(h)','units','normalized'); |
---|
152 | |
---|
153 | subplot(3,3,2) |
---|
154 | for ie = ickd |
---|
155 | plot(-1,-1,cols{ie},'linewidth',1.5); |
---|
156 | hold on |
---|
157 | end |
---|
158 | for ie = ickd |
---|
159 | flux_error = ckd{ie}.flux_up_lw(:,iprof) - ref.flux_up_lw(:,iprof); |
---|
160 | flux_bias = mean(flux_error'); |
---|
161 | flux_ci = std(flux_error').*1.96; |
---|
162 | flux_errmin = min(flux_error,[],2); |
---|
163 | flux_errmax = max(flux_error,[],2); |
---|
164 | pp = mean(p_hl_hPa'); |
---|
165 | h=fill([flux_bias+flux_ci flip(flux_bias-flux_ci)],... |
---|
166 | [pp flip(pp)],'r'); |
---|
167 | set(h,'facecolor',cols{ie}(1),'edgecolor','none','facealpha',alpha(ie)); |
---|
168 | hold on |
---|
169 | plot(flux_bias,pp,cols{ie},'linewidth',1.5); |
---|
170 | end |
---|
171 | set(gca,'ydir','reverse','yscale','log'); |
---|
172 | ylim([0.01 1000]); |
---|
173 | plot([0 0],[0.01 1000],'k:'); |
---|
174 | ylabel('Pressure (hPa)'); |
---|
175 | xlabel('Upwelling flux error (W m^{-2})'); |
---|
176 | xlim(dflux_axis); |
---|
177 | title('Errors in profiles') |
---|
178 | if do_legend |
---|
179 | set(legend(titles{ickd},'location','best'),'fontsize',legfontsize','box','on'); |
---|
180 | end |
---|
181 | text(xloc, yloc, '\bf(b)','units','normalized'); |
---|
182 | |
---|
183 | subplot(3,3,5) |
---|
184 | for ie = ickd |
---|
185 | flux_error = ckd{ie}.flux_dn_lw(:,iprof) - ref.flux_dn_lw(:,iprof); |
---|
186 | flux_bias = mean(flux_error'); |
---|
187 | flux_ci = std(flux_error').*1.96; |
---|
188 | flux_errmin = min(flux_error,[],2); |
---|
189 | flux_errmax = max(flux_error,[],2); |
---|
190 | pp = mean(p_hl_hPa'); |
---|
191 | h=fill([flux_bias+flux_ci flip(flux_bias-flux_ci)],... |
---|
192 | [pp flip(pp)],'r'); |
---|
193 | set(h,'facecolor',cols{ie}(1),'edgecolor','none','facealpha',alpha(ie)); |
---|
194 | hold on |
---|
195 | plot(flux_bias,pp,cols{ie},'linewidth',1.5); |
---|
196 | end |
---|
197 | set(gca,'ydir','reverse','yscale','log'); |
---|
198 | ylim([0.01 1000]); |
---|
199 | plot([0 0],[0.01 1000],'k:'); |
---|
200 | ylabel('Pressure (hPa)'); |
---|
201 | xlabel('Downwelling flux error (W m^{-2})'); |
---|
202 | xlim(dflux_axis); |
---|
203 | text(xloc, yloc, '\bf(e)','units','normalized'); |
---|
204 | |
---|
205 | subplot(3,3,3) |
---|
206 | ileg = 1; |
---|
207 | leg = {}; |
---|
208 | end |
---|
209 | |
---|
210 | for ie = ickd |
---|
211 | mydiff = ckd{ie}.flux_up_lw(1,iprof)-ref.flux_up_lw(1,iprof); |
---|
212 | err = sqrt(mean(mydiff.^2)); |
---|
213 | leg{ileg} = [titles{ie} ' (RMSE=' num2str(err,'%0.2f') ' W m^{-2})']; |
---|
214 | error_save(3,ie) = err; |
---|
215 | error_save(4,ie) = mean(mydiff); |
---|
216 | %plot(-1,-1,cols{ie},'linewidth',1.5); |
---|
217 | ileg = ileg+1; |
---|
218 | if do_plot |
---|
219 | plot(ref.flux_up_lw(1,iprof), mydiff, ... |
---|
220 | [cols{ie}(1) 'o']); |
---|
221 | hold on |
---|
222 | end |
---|
223 | end |
---|
224 | if do_plot |
---|
225 | if do_legend |
---|
226 | set(legend(leg,'location','best'),'fontsize',legfontsize','box','on') |
---|
227 | end |
---|
228 | xlabel('Reference TOA upwelling (W m^{-2})'); |
---|
229 | ylabel('TOA upwelling error (W m^{-2})'); |
---|
230 | xlim([0 400]); |
---|
231 | plot(xlim,[0 0],'k:','linewidth',0.5); |
---|
232 | ylim(dfluxscatter_axis); |
---|
233 | title('Errors at surface and TOA'); |
---|
234 | text(xloc, yloc, '\bf(c)','units','normalized'); |
---|
235 | |
---|
236 | subplot(3,3,6) |
---|
237 | end |
---|
238 | ileg = 1; |
---|
239 | leg = {}; |
---|
240 | for ie = ickd |
---|
241 | mydiff = ckd{ie}.flux_dn_lw(end,iprof)-ref.flux_dn_lw(end,iprof); |
---|
242 | err = sqrt(mean(mydiff.^2)); |
---|
243 | leg{ileg} = [titles{ie} ' (RMSE=' num2str(err,'%0.2f') ' W m^{-2})']; |
---|
244 | error_save(5,ie) = err; |
---|
245 | error_save(6,ie) = mean(mydiff); |
---|
246 | ileg = ileg+1; |
---|
247 | if do_plot |
---|
248 | plot(ref.flux_dn_lw(end,iprof), mydiff, ... |
---|
249 | [cols{ie}(1) 'o']); |
---|
250 | hold on |
---|
251 | end |
---|
252 | end |
---|
253 | if do_plot |
---|
254 | if do_legend |
---|
255 | set(legend(leg,'location','best'),'fontsize',legfontsize') |
---|
256 | end |
---|
257 | xlabel('Reference surface downwelling (W m^{-2})'); |
---|
258 | ylabel('Surface downwelling error (W m^{-2})'); |
---|
259 | plot(xlim,[0 0],'k:','linewidth',0.5); |
---|
260 | ylim(dfluxscatter_axis); |
---|
261 | text(xloc, yloc, '\bf(f)','units','normalized'); |
---|
262 | |
---|
263 | h=subplot(3,3,9); |
---|
264 | set(h,'visible','off') |
---|
265 | |
---|
266 | xstart = -0.15; |
---|
267 | text(xstart,0.95,['\bfScenario: ' scenario_title],'fontsize',12); |
---|
268 | |
---|
269 | if ickd_stats > 0 |
---|
270 | for ie = ickd_stats |
---|
271 | text(xstart,0.85,['\bfCKD model: ' titles{ie}],'fontsize',12); |
---|
272 | ystart = 0.75; %-0.4.*(ie-1); |
---|
273 | %text(xstart,ystart,['\bfRMS errors in ' titles{ie} ':'],'units','normalized'); |
---|
274 | text(xstart,ystart,['Bias TOA upwelling: ' num2str(error_save(4,ie),'%0.2f') ' W m^{-2}'],'units','normalized'); |
---|
275 | text(xstart,ystart-0.1,['Bias surface downwelling: ' num2str(error_save(6,ie),'%0.2f') ' W m^{-2}'],'units','normalized'); |
---|
276 | text(xstart,ystart-0.2,['RMSE TOA upwelling: ' num2str(error_save(3,ie),'%0.2f') ' W m^{-2}'],'units','normalized'); |
---|
277 | text(xstart,ystart-0.3,['RMSE surface downwelling: ' num2str(error_save(5,ie),'%0.2f') ' W m^{-2}'],'units','normalized'); |
---|
278 | text(xstart,ystart-0.4,['RMSE heating rate (0.02-4 hPa): ' num2str(error_save(1,ie),'%0.3f') ' K d^{-1}'],'units','normalized'); |
---|
279 | text(xstart,ystart-0.5,['RMSE heating rate (4-1100 hPa): ' num2str(error_save(2,ie),'%0.3f') ' K d^{-1}'],'units','normalized'); |
---|
280 | |
---|
281 | end |
---|
282 | end |
---|
283 | |
---|
284 | drawnow |
---|
285 | end |
---|
286 | %print('-dpng','-painters','-r100',[combined_titles 'evaluation1_fluxes_' scenario '.png']); |
---|
287 | %end |
---|
288 | |
---|
289 | if nargin > 0 |
---|
290 | stats.toa_up_bias = error_save(4,:); |
---|
291 | stats.surf_dn_bias= error_save(6,:); |
---|
292 | stats.toa_up_rmse = error_save(3,:); |
---|
293 | stats.surf_dn_rmse= error_save(5,:); |
---|
294 | stats.heating_rate_low_rmse=error_save(1,:); |
---|
295 | stats.heating_rate_high_rmse=error_save(2,:); |
---|
296 | end |
---|
297 | |
---|