1 ;********************************************************
2 ;using model biome vlass
4 ; required command line input parameters:
5 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
7 ; histogram normalized by rain and compute correleration
8 ;**************************************************************
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
12 ;**************************************************************
13 procedure set_line(lines:string,nline:integer,newlines:string)
15 ; add line to ascci/html file
17 nnewlines = dimsizes(newlines)
18 if(nline+nnewlines-1.ge.dimsizes(lines))
19 print("set_line: bad index, not setting anything.")
22 lines(nline:nline+nnewlines-1) = newlines
23 ; print ("lines = " + lines(nline:nline+nnewlines-1))
24 nline = nline + nnewlines
27 ;**************************************************************
34 ;---------------------------------------------------
40 model_name1 = "i01.06cn"
41 model_name2 = "i01.10cn"
43 ;---------------------------------------------------
44 ; get biome data: model
46 biome_name_mod = "Model PFT Class"
48 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
49 film = "class_pft_"+model_grid+".nc"
50 fm = addfile(dirm+film,"r")
52 classmod = fm->CLASS_PFT
56 ; model data has 17 land-type classes
60 ;--------------------------------------------------
61 ; get model data: landmask, landfrac and area
63 dirm = "/fis/cgd/cseg/people/jeff/surface_data/"
65 fm = addfile (dirm+film,"r")
67 landmask = fm->landmask
68 landfrac = fm->landfrac
73 ; change area from km**2 to m**2
76 ;----------------------------------------------------
77 ; read data: time series, model
79 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
80 film = model_name2 + "_Fire_C_1979-2004_monthly.nc"
81 fm = addfile (dirm+film,"r")
83 data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)
87 ; Units for these variables are:
90 ; change unit to g C/m^2/month
92 nsec_per_month = 60*60*24*30
94 data_mod = data_mod * nsec_per_month
96 data_mod@unit = "gC/m2/month"
97 ;----------------------------------------------------
98 ; read data: time series, observed
100 dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
101 film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"
102 fm = addfile (dirm+film,"r")
104 data_ob = fm->FIRE_C(0:7,:,:,:)
110 ; Units for these variables are:
113 ;---------------------------------------------------
114 ; take into account landfrac
116 area = area * landfrac
117 data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))
118 data_ob = data_ob * conform(data_ob, landfrac, (/2,3/))
122 ;-------------------------------------------------------------
125 ; column (not including header column)
127 col_head = (/"Observed Fire_Flux (PgC/yr)" \
128 ,"Model Fire_Flux (PgC/yr)" \
129 ,"Correlation Coefficient" \
130 ,"Ratio model/observed" \
135 ncol = dimsizes(col_head)
137 ; row (not including header row)
139 ; using model biome class:
140 row_head = (/"Not Vegetated" \
141 ,"Needleleaf Evergreen Temperate Tree" \
142 ,"Needleleaf Evergreen Boreal Tree" \
143 ; ,"Needleleaf Deciduous Boreal Tree" \
144 ,"Broadleaf Evergreen Tropical Tree" \
145 ,"Broadleaf Evergreen Temperate Tree" \
146 ,"Broadleaf Deciduous Tropical Tree" \
147 ,"Broadleaf Deciduous Temperate Tree" \
148 ; ,"Broadleaf Deciduous Boreal Tree" \
149 ; ,"Broadleaf Evergreen Shrub" \
150 ,"Broadleaf Deciduous Temperate Shrub" \
151 ,"Broadleaf Deciduous Boreal Shrub" \
153 ,"C3 Non-Arctic Grass" \
159 nrow = dimsizes(row_head)
161 ; arrays to be passed to table.
162 text = new ((/nrow, ncol/),string )
164 ;*****************************************************************
166 ;*****************************************************************
168 x = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))
169 data_mod_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
172 x = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))
173 data_ob_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
176 ;----------------------------------------------------
177 ; compute correlation coef
179 landmask_1d = ndtooned(landmask)
180 data_mod_1d = ndtooned(data_mod_m)
181 data_ob_1d = ndtooned(data_ob_m)
183 good = ind(landmask_1d .gt. 0.)
185 cc = esccr(data_mod_1d(good),data_ob_1d(good),0)
192 ;----------------------------------------------------
197 Mscore1 = cc * cc * score_max
199 M_global = sprintf("%.2f", Mscore1)
201 ;----------------------------------------------------
204 resg = True ; Use plot options
205 resg@cnFillOn = True ; Turn on color fill
206 resg@gsnSpreadColors = True ; use full colormap
207 resg@cnLinesOn = False ; Turn off contourn lines
208 resg@mpFillOn = False ; Turn off map fill
209 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
211 ;----------------------------------------------------
212 ; global contour: model vs ob
214 plot_name = "global_model_vs_ob"
216 wks = gsn_open_wks (plot_type,plot_name)
217 gsn_define_colormap(wks,"gui_default")
219 plot=new(3,graphic) ; create graphic array
221 resg@gsnFrame = False ; Do not draw plot
222 resg@gsnDraw = False ; Do not advance frame
224 ;----------------------
225 ; plot correlation coef
228 gRes@txFontHeightF = 0.02
231 correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"
233 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
235 ;-----------------------
238 data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)
241 resg@tiMainString = title
243 resg@cnMinLevelValF = 1.
244 resg@cnMaxLevelValF = 10.
245 resg@cnLevelSpacingF = 1.
247 plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg)
249 ;-----------------------
252 data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)
254 title = "Model "+ model_name
255 resg@tiMainString = title
257 resg@cnMinLevelValF = 1.
258 resg@cnMaxLevelValF = 10.
259 resg@cnLevelSpacingF = 1.
261 plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg)
263 ;-----------------------
266 resg@cnMinLevelValF = -8.
267 resg@cnMaxLevelValF = 2.
268 resg@cnLevelSpacingF = 1.
271 zz = data_mod_m - data_ob_m
272 title = "Model_"+model_name+" - Observed"
273 resg@tiMainString = title
275 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
279 pres = True ; panel plot mods desired
280 pres@gsnMaximize = True ; fill the page
282 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
284 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
285 "rm "+plot_name+"."+plot_type)
294 resg@gsnFrame = True ; Do advance frame
295 resg@gsnDraw = True ; Do draw plot
297 ;*******************************************************************
298 ; (B) Time series : per biome
299 ;*******************************************************************
303 dsizes = dimsizes(data_mod)
306 ntime = nyear * nmonth
311 ;-------------------------------------------
312 ; Calculate "nice" bins for binning the data
314 ; using model biome class
317 range = fspan(0,nclass,nclass+1)
319 ; Use this range information to grab all the values in a
320 ; particular range, and then take an average.
322 nx = dimsizes(range) - 1
324 ;-------------------------------------------
327 ; using observed biome class
328 ; base = ndtooned(classob)
329 ; using model biome class
330 base = ndtooned(classmod)
334 area_bin = new((/nx/),float)
335 yvalues = new((/ntime,data_n,nx/),float)
337 ; Loop through each range, using base.
341 if (i.ne.(nx-1)) then
342 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
344 idx = ind(base.ge.range(i))
346 ;---------------------
349 data = ndtooned(area)
351 if (.not.any(ismissing(idx))) then
352 area_bin(i) = sum(data(idx))
354 area_bin(i) = area_bin@_FillValue
357 ;#############################################################
358 ; using model biome class:
359 ; set the following 4 classes to _FillValue:
360 ; (3)Needleleaf Deciduous Boreal Tree,
361 ; (8)Broadleaf Deciduous Boreal Tree,
362 ; (9)Broadleaf Evergreen Shrub,
365 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
366 area_bin(i) = area_bin@_FillValue
368 ;#############################################################
372 ;---------------------
373 ; for data_mod and data_ob
384 data = ndtooned(data_ob(m,k,:,:))
388 data = ndtooned(data_mod(m,k,:,:))
393 if (.not.any(ismissing(idx))) then
394 yvalues(t,n,i) = avg(data(idx))
396 yvalues(t,n,i) = yvalues@_FillValue
399 ;#############################################################
400 ; using model biome class:
401 ; set the following 4 classes to _FillValue:
402 ; (3)Needleleaf Deciduous Boreal Tree,
403 ; (8)Broadleaf Deciduous Boreal Tree,
404 ; (9)Broadleaf Evergreen Shrub,
407 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
408 yvalues(t,n,i) = yvalues@_FillValue
410 ;#############################################################
425 ;----------------------------------------------------------------
428 good = ind(.not.ismissing(area_bin))
430 area_g = area_bin(good)
432 n_biome = dimsizes(good)
434 ;----------------------------------------------------------------
435 ; data for tseries plot
437 yvalues_g = new((/ntime,data_n,n_biome/),float)
439 yvalues_g@units = "TgC/month"
441 ; change unit to Tg C/month
442 ; change unit from g to Tg (Tera gram)
445 yvalues_g = yvalues(:,:,good) * conform(yvalues_g,area_g,2) * factor_unit
449 ;-------------------------------------------------------------------
450 ; general settings for line plot
453 res@xyDashPatterns = (/0,0/) ; make lines solid
454 res@xyLineThicknesses = (/2.0,2.0/) ; make lines thicker
455 res@xyLineColors = (/"blue","red"/) ; line color
457 res@trXMinF = year_start
458 res@trXMaxF = year_end + 1
460 res@vpHeightF = 0.4 ; change aspect ratio of plot
464 res@tiMainFontHeightF = 0.025 ; size of title
466 res@tmXBFormat = "f" ; not to add trailing zeros
468 ; res@gsnMaximize = True
470 ;----------------------------------------------
471 ; Add a boxed legend using the simple method
473 res@pmLegendDisplayMode = "Always"
474 ; res@pmLegendWidthF = 0.1
475 res@pmLegendWidthF = 0.08
476 res@pmLegendHeightF = 0.06
477 res@pmLegendOrthogonalPosF = -1.17
478 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
479 ; res@pmLegendOrthogonalPosF = -0.30 ;(downward)
481 ; res@pmLegendParallelPosF = 0.18
482 res@pmLegendParallelPosF = 0.23 ;(rightward)
483 res@pmLegendParallelPosF = 0.73 ;(rightward)
484 res@pmLegendParallelPosF = 0.83 ;(rightward)
486 ; res@lgPerimOn = False
487 res@lgLabelFontHeightF = 0.015
488 res@xyExplicitLegendLabels = (/"observed",model_name/)
490 ;*******************************************************************
491 ; (A) time series plot: monthly ( 2 lines per plot)
492 ;*******************************************************************
494 ; x-axis in time series plot
496 timeI = new((/ntime/),integer)
497 timeF = new((/ntime/),float)
498 timeI = ispan(1,ntime,1)
499 timeF = year_start + (timeI-1)/12.
500 timeF@long_name = "year"
502 plot_data = new((/2,ntime/),float)
503 plot_data@long_name = "TgC/month"
505 ;----------------------------------------------
506 ; time series plot : per biome
510 plot_name = "monthly_biome_"+ m
512 wks = gsn_open_wks (plot_type,plot_name)
514 title = "Fire : "+ row_head(m)
515 res@tiMainString = title
517 plot_data(0,:) = yvalues_g(:,0,m)
518 plot_data(1,:) = yvalues_g(:,1,m)
520 plot = gsn_csm_xy(wks,timeF,plot_data,res)
522 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
523 "rm "+plot_name+"."+plot_type)
530 ;------------------------------------------
531 ; data for table : per biome
533 ; unit change from TgC/month to PgC/month
538 tmp_ob = new((/ntime/),float)
539 tmp_mod = new((/ntime/),float)
541 total_ob = new((/n_biome/),float)
542 total_mod = new((/n_biome/),float)
543 Mscore2 = new((/n_biome/),float)
547 tmp_ob = yvalues_g(:,0,m)
548 tmp_mod = yvalues_g(:,1,m)
550 total_ob(m) = avg(month_to_annual(tmp_ob, 0)) * unit_factor
551 total_mod(m) = avg(month_to_annual(tmp_mod,0)) * unit_factor
553 cc = esccr(tmp_mod,tmp_ob,0)
555 ratio = total_mod(m)/total_ob(m)
557 good = ind(tmp_ob .ne. 0. .and. tmp_mod .ne. 0.)
559 bias = sum( abs( tmp_mod(good)-tmp_ob(good) )/( abs(tmp_mod(good))+abs(tmp_ob(good)) ) )
560 Mscore2(m) = (1.- (bias/dimsizes(good)))*score_max
564 text(m,0) = sprintf("%.2f",total_ob(m))
565 text(m,1) = sprintf("%.2f",total_mod(m))
566 text(m,2) = sprintf("%.2f",cc)
567 text(m,3) = sprintf("%.2f",ratio)
568 text(m,4) = sprintf("%.2f",Mscore2(m))
569 text(m,5) = "<a href=./monthly_biome_"+m+".png>model_vs_ob</a>"
575 ;--------------------------------------------
576 ; time series plot: all biome
578 plot_name = "monthly_global"
580 wks = gsn_open_wks (plot_type,plot_name)
582 title = "Fire : "+ row_head(n_biome)
583 res@tiMainString = title
586 plot_data(0,k) = sum(yvalues_g(k,0,:))
587 plot_data(1,k) = sum(yvalues_g(k,1,:))
590 plot = gsn_csm_xy(wks,timeF,plot_data,res)
592 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
593 "rm "+plot_name+"."+plot_type)
598 ;------------------------------------------
599 ; data for table : global
601 tmp_ob = ndtooned(yvalues_g(:,0,:))
602 tmp_mod = ndtooned(yvalues_g(:,1,:))
604 cc = esccr(tmp_mod,tmp_ob,0)
606 ratio = sum(total_mod)/sum(total_ob)
608 text(nrow-1,0) = sprintf("%.2f",sum(total_ob))
609 text(nrow-1,1) = sprintf("%.2f",sum(total_mod))
610 text(nrow-1,2) = sprintf("%.2f",cc)
611 text(nrow-1,3) = sprintf("%.2f",ratio)
612 text(nrow-1,4) = sprintf("%.2f",avg(Mscore2))
613 text(nrow-1,5) = "<a href=./monthly_global.png>model_vs_ob</a>"
615 ;**************************************************
617 ;**************************************************
619 header_text = "<H1>Fire Emission (1997-2004): Model "+model_name+"</H1>"
621 header = (/"<HTML>" \
623 ,"<TITLE>CLAMP metrics</TITLE>" \
630 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
632 ," <th bgcolor=DDDDDD >Biome Type</th>" \
633 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
634 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
635 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
636 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
637 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
638 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
641 table_footer = "</table>"
645 lines = new(50000,string)
648 set_line(lines,nline,header)
649 set_line(lines,nline,table_header)
650 ;-----------------------------------------------
654 set_line(lines,nline,row_header)
664 set_line(lines,nline,"<th>"+txt0+"</th>")
665 set_line(lines,nline,"<th>"+txt1+"</th>")
666 set_line(lines,nline,"<th>"+txt2+"</th>")
667 set_line(lines,nline,"<th>"+txt3+"</th>")
668 set_line(lines,nline,"<th>"+txt4+"</th>")
669 set_line(lines,nline,"<th>"+txt5+"</th>")
670 set_line(lines,nline,"<th>"+txt6+"</th>")
672 set_line(lines,nline,row_footer)
674 ;-----------------------------------------------
675 set_line(lines,nline,table_footer)
676 set_line(lines,nline,footer)
678 ; Now write to an HTML file.
680 output_html = "table_fire.html"
682 idx = ind(.not.ismissing(lines))
683 if(.not.any(ismissing(idx))) then
684 asciiwrite(output_html,lines(idx))