source: LMDZ6/branches/contrails/libf/phylmd/ecrad/test/ckdmip/evaluate_ckd_sw_fluxes.m

Last change on this file was 4773, checked in by idelkadi, 12 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


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