Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;**************************************************************
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
5 ;**************************************************************
6 procedure set_line(lines:string,nline:integer,newlines:string)
8 ; add line to ascci/html file
10 nnewlines = dimsizes(newlines)
11 if(nline+nnewlines-1.ge.dimsizes(lines))
12 print("set_line: bad index, not setting anything.")
15 lines(nline:nline+nnewlines-1) = newlines
16 ; print ("lines = " + lines(nline:nline+nnewlines-1))
17 nline = nline + nnewlines
20 ;**************************************************************
27 ;----------------------------------------------------------
28 ; edit current model for movel1_vs_model2
30 if (isvar("compare")) then
31 html_name2 = compare+"/table.html"
32 html_new2 = html_name2 +".new"
35 ;----------------------------------------------------------
36 ; edit table.html for current model
38 html_name = model_name+"/table.html"
39 html_new = html_name +".new"
41 ;----------------------------------------------------------
42 ; get biome data: model
44 biome_name_mod = "Model PFT Class"
46 film_c = "class_pft_"+ model_grid +".nc"
47 fm_c = addfile (dirs+film_c,"r")
48 classmod = fm_c->CLASS_PFT
52 ; model data has 17 land-type classes
55 ;------------------------------------
56 ; get model landfrac and area
58 film_l = "lnd_"+ model_grid +".nc"
59 fm_l = addfile (dirs+film_l,"r")
60 landfrac = fm_l->landfrac
65 ; change area from km**2 to m**2
68 ;-------------------------------------
69 ; take into account landfrac
71 ; area = area * landfrac
75 ;---------------------------------------------------
76 ; read data: model, group 1
78 fm = addfile (dirm+film4,"r")
86 VegC = leafc + woodc + frootc
91 LiCwC = litterc + cwdc
96 ;---------------------------------------------------
97 ; read data: model, group 2
99 fm = addfile (dirm+film5,"r")
106 ;---------------------------------------------------
107 ; Units for these variables are:
118 nsec_per_year = 60*60*24*365
120 ; change unit to g C/m^2/year
122 NPP1 = NPP1 * nsec_per_year * conform(NPP1,landfrac,(/1,2/))
123 NPP2 = NPP2 * nsec_per_year * conform(NPP2,landfrac,(/1,2/))
124 NEE2 = NEE2 * nsec_per_year * conform(NEE2,landfrac,(/1,2/))
125 GPP2 = GPP2 * nsec_per_year * conform(GPP2,landfrac,(/1,2/))
127 VegC = VegC * conform(VegC,landfrac,(/1,2/))
128 LiCwC = LiCwC * conform(LiCwC,landfrac,(/1,2/))
129 SoilC = SoilC * conform(SoilC,landfrac,(/1,2/))
133 ;*******************************************************************
134 ; Calculate "nice" bins for binning the data in equally spaced ranges
135 ;********************************************************************
137 ; using model biome class
140 range = fspan(0,nclass,nclass+1)
143 ; Use this range information to grab all the values in a
144 ; particular range, and then take an average.
146 nx = dimsizes(range) - 1
148 ;==============================
150 ;==============================
152 ; using observed biome class
153 ; base = ndtooned(classob)
154 ; using model biome class
155 base = ndtooned(classmod)
157 area_1d = ndtooned(area)
161 yvalues = new((/data_n,nx/),float) ; (per m2)
162 yvalues_t = new((/data_n,nx/),float) ; (per biome)
164 ; Loop through each range, using base.
168 if (i.ne.(nx-1)) then
169 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
171 idx = ind(base.ge.range(i))
177 data = ndtooned(area)
181 data = ndtooned(NPP1)
185 data = ndtooned(VegC)
189 data = ndtooned(LiCwC)
193 data = ndtooned(SoilC)
197 data = ndtooned(NPP2)
201 data = ndtooned(NEE2)
205 data = ndtooned(GPP2)
208 ; Calculate sum and average
210 if (.not.any(ismissing(idx))) then
212 yvalues(n,i) = sum(data(idx))
213 yvalues_t(n,i) = sum(data(idx))
215 yvalues(n,i) = avg(data(idx))
216 yvalues_t(n,i) = sum(data(idx)*area_1d(idx))
219 yvalues(n,i) = yvalues@_FillValue
220 yvalues_t(n,i) = yvalues@_FillValue
223 ;#############################################################
224 ; using model biome class:
225 ; set the following 4 classes to _FillValue:
226 ; (3)Needleleaf Deciduous Boreal Tree,
227 ; (8)Broadleaf Deciduous Boreal Tree,
228 ; (9)Broadleaf Evergreen Shrub,
231 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
232 yvalues(n,i) = yvalues@_FillValue
233 yvalues_t(n,i) = yvalues@_FillValue
235 ;#############################################################
253 ;----------------------------------------------------------------
256 good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
259 area_g = yvalues(0,good)
260 NPP1_g = yvalues(1,good)
261 VegC_g = yvalues(2,good)
262 LiCwC_g = yvalues(3,good)
263 SoilC_g = yvalues(4,good)
264 NPP2_g = yvalues(5,good)
265 NEE2_g = yvalues(6,good)
266 GPP2_g = yvalues(7,good)
268 NPP_ratio = NPP2_g/NPP1_g
270 n_biome = dimsizes(NPP1_g)
272 ;-----------------------------------------------------------------
275 ; change unit from g to Pg (Peta gram)
278 NPP1_t = yvalues_t(1,good) * factor_unit
279 VegC_t = yvalues_t(2,good) * factor_unit
280 LiCwC_t = yvalues_t(3,good) * factor_unit
281 SoilC_t = yvalues_t(4,good) * factor_unit
282 NEE2_t = yvalues_t(6,good) * factor_unit
283 GPP2_t = yvalues_t(7,good) * factor_unit
288 ;-------------------------------------------------------------
291 ; column (not including header column)
293 col_head = (/"Area (1.e12m2)" \
296 ,"Litter+CWD (gC/m2)" \
303 ncol = dimsizes(col_head)
305 ; row (not including header row)
307 ; using model biome class:
308 row_head = (/"Not Vegetated" \
309 ,"Needleleaf Evergreen Temperate Tree" \
310 ,"Needleleaf Evergreen Boreal Tree" \
311 ; ,"Needleleaf Deciduous Boreal Tree" \
312 ,"Broadleaf Evergreen Tropical Tree" \
313 ,"Broadleaf Evergreen Temperate Tree" \
314 ,"Broadleaf Deciduous Tropical Tree" \
315 ,"Broadleaf Deciduous Temperate Tree" \
316 ; ,"Broadleaf Deciduous Boreal Tree" \
317 ; ,"Broadleaf Evergreen Shrub" \
318 ,"Broadleaf Deciduous Temperate Shrub" \
319 ,"Broadleaf Deciduous Boreal Shrub" \
321 ,"C3 Non-Arctic Grass" \
327 nrow = dimsizes(row_head)
329 ; arrays to be passed to table.
330 text = new ((/nrow, ncol/),string )
333 text(i,0) = sprintf("%.1f",area_g(i)*1.e-12)
334 text(i,1) = sprintf("%.1f",NPP1_g(i))
335 text(i,2) = sprintf("%.1f",VegC_g(i))
336 text(i,3) = sprintf("%.1f",LiCwC_g(i))
337 text(i,4) = sprintf("%.1f",SoilC_g(i))
338 text(i,5) = sprintf("%.2f",NPP_ratio(i))
339 text(i,6) = sprintf("%.1f",NEE2_g(i))
340 text(i,7) = sprintf("%.1f",GPP2_g(i))
343 ;-------------------------------------------------------
346 header_text = "<H1>NEE and Carbon Stocks and Fluxes: Model "+model_name+"</H1>"
348 header = (/"<HTML>" \
350 ,"<TITLE>CLAMP metrics</TITLE>" \
357 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
359 ," <th bgcolor=DDDDDD >Biome Type</th>" \
360 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
361 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
362 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
363 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
364 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
365 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
366 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
367 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
370 table_footer = "</table>"
374 lines = new(50000,string)
377 set_line(lines,nline,header)
378 set_line(lines,nline,table_header)
380 ;----------------------------
384 set_line(lines,nline,row_header)
396 set_line(lines,nline,"<th>"+txt0+"</th>")
397 set_line(lines,nline,"<th>"+txt1+"</th>")
398 set_line(lines,nline,"<th>"+txt2+"</th>")
399 set_line(lines,nline,"<th>"+txt3+"</th>")
400 set_line(lines,nline,"<th>"+txt4+"</th>")
401 set_line(lines,nline,"<th>"+txt5+"</th>")
402 set_line(lines,nline,"<th>"+txt6+"</th>")
403 set_line(lines,nline,"<th>"+txt7+"</th>")
404 set_line(lines,nline,"<th>"+txt8+"</th>")
406 set_line(lines,nline,row_footer)
408 ;----------------------------
409 set_line(lines,nline,table_footer)
410 set_line(lines,nline,footer)
412 ; Now write to an HTML file.
414 output_html = "table_per_m2.html"
416 idx = ind(.not.ismissing(lines))
417 if(.not.any(ismissing(idx))) then
418 asciiwrite(output_html,lines(idx))
428 delete (table_header)
430 ;-----------------------------------------------------------------
433 ; column (not including header column)
435 col_head = (/"NPP (PgC/yr)" \
437 ,"Litter+CWD (PgC)" \
446 ncol = dimsizes(col_head)
448 ; row (not including header row)
450 ; using model biome class:
451 row_head = (/"Not Vegetated" \
452 ,"Needleleaf Evergreen Temperate Tree" \
453 ,"Needleleaf Evergreen Boreal Tree" \
454 ; ,"Needleleaf Deciduous Boreal Tree" \
455 ,"Broadleaf Evergreen Tropical Tree" \
456 ,"Broadleaf Evergreen Temperate Tree" \
457 ,"Broadleaf Deciduous Tropical Tree" \
458 ,"Broadleaf Deciduous Temperate Tree" \
459 ; ,"Broadleaf Deciduous Boreal Tree" \
460 ; ,"Broadleaf Evergreen Shrub" \
461 ,"Broadleaf Deciduous Temperate Shrub" \
462 ,"Broadleaf Deciduous Boreal Shrub" \
464 ,"C3 Non-Arctic Grass" \
470 nrow = dimsizes(row_head)
472 ; arrays to be passed to table.
473 text = new ((/nrow, ncol/),string )
476 text(i,0) = sprintf("%.1f",NPP1_t(i))
477 text(i,1) = sprintf("%.1f",VegC_t(i))
478 text(i,2) = sprintf("%.1f",LiCwC_t(i))
479 text(i,3) = sprintf("%.1f",SoilC_t(i))
480 text(i,4) = sprintf("%.1f",NEE2_t(i))
481 text(i,5) = sprintf("%.1f",GPP2_t(i))
482 text(i,6) = "<a href=./NPP_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NPP_annual_biome_"+i+".png>annual_plot</a>"
483 text(i,7) = "<a href=./NEE_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NEE_annual_biome_"+i+".png>annual_plot</a>"
487 text(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))
488 text(nrow-1,1) = sprintf("%.1f",sum(VegC_t))
489 text(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))
490 text(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))
491 text(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))
492 text(nrow-1,5) = sprintf("%.1f",sum(GPP2_t))
493 text(nrow-1,6) = "<a href=./NPP_monthly_global.png>monthly_plot</a> <br> <a href=./NPP_annual_global.png>annual_plot</a>"
494 text(nrow-1,7) = "<a href=./NEE_monthly_global.png>monthly_plot</a> <br> <a href=./NEE_annual_global.png>annual_plot</a>"
495 text(nrow-1,8) = "--"
497 ;**************************************************
499 ;**************************************************
501 header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome): Model "+model_name+"</H1>"
503 header = (/"<HTML>" \
505 ,"<TITLE>CLAMP metrics</TITLE>" \
512 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
514 ," <th bgcolor=DDDDDD >Biome Type</th>" \
515 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
516 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
517 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
518 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
519 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
520 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
521 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
522 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
523 ," <th bgcolor=DDDDDD >"+col_head(8)+"</th>" \
526 table_footer = "</table>"
530 lines = new(50000,string)
533 set_line(lines,nline,header)
534 set_line(lines,nline,table_header)
535 ;-----------------------------------------------
539 set_line(lines,nline,row_header)
552 set_line(lines,nline,"<th>"+txt0+"</th>")
553 set_line(lines,nline,"<th>"+txt1+"</th>")
554 set_line(lines,nline,"<th>"+txt2+"</th>")
555 set_line(lines,nline,"<th>"+txt3+"</th>")
556 set_line(lines,nline,"<th>"+txt4+"</th>")
557 set_line(lines,nline,"<th>"+txt5+"</th>")
558 set_line(lines,nline,"<th>"+txt6+"</th>")
559 set_line(lines,nline,"<th>"+txt7+"</th>")
560 set_line(lines,nline,"<th>"+txt8+"</th>")
561 set_line(lines,nline,"<th>"+txt9+"</th>")
563 set_line(lines,nline,row_footer)
565 ;-----------------------------------------------
566 set_line(lines,nline,table_footer)
567 set_line(lines,nline,footer)
569 ; Now write to an HTML file.
571 output_html = "table_per_biome.html"
573 idx = ind(.not.ismissing(lines))
574 if(.not.any(ismissing(idx))) then
575 asciiwrite(output_html,lines(idx))
582 ;---------------------------------------------------
583 ; read model data, time series:
585 fm = addfile (dirm+film7,"r")
589 ;Fire = fm->COL_FIRE_CLOSS
593 ; Units for these variables are:
599 nsec_per_month = 60*60*24*30
601 ; change unit to g C/m^2/month
603 NPP3 = NPP3 * nsec_per_month * conform(NPP3,landfrac,(/2,3/))
604 NEE3 = NEE3 * nsec_per_month * conform(NEE3,landfrac,(/2,3/))
605 ;Fire = Fire * nsec_per_month * conform(Fire,landfrac,(/2,3/))
610 dsizes = dimsizes(NPP3)
613 ntime = nyear * nmonth
618 ;*******************************************************************
619 ; Calculate "nice" bins for binning the data in equally spaced ranges
620 ;********************************************************************
622 ; using model biome class
625 range = fspan(0,nclass,nclass+1)
628 ; Use this range information to grab all the values in a
629 ; particular range, and then take an average.
631 nx = dimsizes(range) - 1
633 ;==============================
635 ;==============================
637 ; using observed biome class
638 ; base = ndtooned(classob)
639 ; using model biome class
640 base = ndtooned(classmod)
644 yvalues = new((/ntime,data_n,nx/),float)
646 ; Loop through each range, using base.
650 if (i.ne.(nx-1)) then
651 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
653 idx = ind(base.ge.range(i))
665 data = ndtooned(NPP3(m,k,:,:))
669 data = ndtooned(NEE3(m,k,:,:))
673 ; data = ndtooned(Fire(m,k,:,:))
678 if (.not.any(ismissing(idx))) then
679 yvalues(t,n,i) = sum(data(idx)*area_1d(idx))
681 yvalues(t,n,i) = yvalues@_FillValue
684 ;#############################################################
685 ; using model biome class:
686 ; set the following 4 classes to _FillValue:
687 ; (3)Needleleaf Deciduous Boreal Tree,
688 ; (8)Broadleaf Deciduous Boreal Tree,
689 ; (9)Broadleaf Evergreen Shrub,
692 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
693 yvalues(t,n,i) = yvalues@_FillValue
695 ;#############################################################
711 ;----------------------------------------------------------------
712 ; data for tseries plot
714 yvalues_g = new((/ntime,data_n,n_biome/),float)
716 yvalues_g@units = "TgC/month"
718 ; change unit to Tg C/month
719 ; change unit from g to Tg (Tera gram)
722 yvalues_g = yvalues(:,:,good) * factor_unit
724 ;*******************************************************************
725 ; general settings for line plot
726 ;*******************************************************************
730 res@xyDashPatterns = (/0/) ; make lines solid
731 res@xyLineThicknesses = (/2.0/) ; make lines thicker
732 res@xyLineColors = (/"blue"/) ; line color
734 res@trXMinF = year_start
735 res@trXMaxF = year_end + 1
737 res@vpHeightF = 0.4 ; change aspect ratio of plot
741 ; res@gsnMaximize = True
743 ;*******************************************************************
744 ; (A) 1 component in each biome: monthly
745 ;*******************************************************************
747 ; component = (/"NPP","NEE","Fire"/)
748 component = (/"NPP","NEE"/)
750 ; for x-axis in xyplot
752 timeI = new((/ntime/),integer)
753 timeF = new((/ntime/),float)
754 timeI = ispan(1,ntime,1)
755 timeF = year_start + (timeI-1)/12.
756 timeF@long_name = "year"
758 plot_data = new((/ntime/),float)
759 plot_data@long_name = "TgC/month"
764 plot_name = component(n)+"_monthly_biome_"+ m
766 wks = gsn_open_wks (plot_type,plot_name)
768 title = component(n)+ ": "+ row_head(m)
769 res@tiMainString = title
770 res@tiMainFontHeightF = 0.025
772 plot_data(:) = yvalues_g(:,n,m)
774 plot=gsn_csm_xy(wks,timeF,plot_data,res)
779 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
780 "rm "+plot_name+"."+plot_type)
786 plot_name = component(n)+"_monthly_global"
788 wks = gsn_open_wks (plot_type,plot_name)
790 title = component(n)+ ": Global"
791 res@tiMainString = title
792 res@tiMainFontHeightF = 0.025
795 plot_data(k) = sum(yvalues_g(k,n,:))
798 plot=gsn_csm_xy(wks,timeF,plot_data,res)
803 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
804 "rm "+plot_name+"."+plot_type)
811 ;*******************************************************************
812 ; (B) 1 component in each biome: annually
813 ;*******************************************************************
815 yvalues_a = new((/nyear,data_n,n_biome/),float)
818 yvalues_g!2 = "record"
820 yvalues_a = month_to_annual(yvalues_g,0)
824 ; for x-axis in xyplot
826 timeI = new((/nyear/),integer)
827 timeF = new((/nyear/),float)
828 timeI = ispan(1,nyear,1)
829 timeF = year_start + (timeI-1)
830 timeF@long_name = "year"
832 plot_data = new((/nyear/),float)
833 plot_data@long_name = "TgC/year"
838 plot_name = component(n)+"_annual_biome_"+ m
840 wks = gsn_open_wks (plot_type,plot_name)
842 title = component(n)+ ": "+ row_head(m)
843 res@tiMainString = title
844 res@tiMainFontHeightF = 0.025
846 plot_data(:) = yvalues_a(:,n,m)
848 plot=gsn_csm_xy(wks,timeF,plot_data,res)
853 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
854 "rm "+plot_name+"."+plot_type)
860 plot_name = component(n)+"_annual_global"
862 wks = gsn_open_wks (plot_type,plot_name)
864 title = component(n)+ ": Global"
865 res@tiMainString = title
866 res@tiMainFontHeightF = 0.025
869 plot_data(k) = sum(yvalues_a(k,n,:))
872 plot=gsn_csm_xy(wks,timeF,plot_data,res)
877 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
878 "rm "+plot_name+"."+plot_type)
881 ;****************************************
882 ; output plot and html
883 ;****************************************
884 output_dir = model_name+"/carbon_sink"
886 system("mv *.png *.html " + output_dir)
887 ;****************************************