Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
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 ;************************************************
36 ;************************************************
41 model_name1 = "i01.06cn"
42 model_name2 = "i01.10cn"
44 ;---------------------------------------------------
45 ; get biome data: model
47 biome_name_mod = "Model PFT Class"
49 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
50 film = "class_pft_"+model_grid+".nc"
52 fm = addfile(dirm+film,"r")
54 classmod = fm->CLASS_PFT
58 ; model data has 17 land-type classes
62 ;--------------------------------------------------
63 ; get model data: landfrac and area
65 dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/"
67 fm_l = addfile (dirm_l+film_l,"r")
69 landfrac = fm_l->landfrac
72 ; change area from km**2 to m**2
75 ;---------------------------------------------------
76 ; take into account landfrac
78 area = area * landfrac
82 ;---------------------------------------------------
83 ; read data: model, group 1
85 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
86 film = model_name1 + "_1980-2004_ANN_climo.nc"
87 fm = addfile (dirm+film,"r")
95 VegC = leafc + woodc + frootc
100 LiCwC = litterc + cwdc
105 ;---------------------------------------------------
106 ; read data: model, group 2
108 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
109 film = model_name2 + "_1990-2004_ANN_climo.nc"
110 fm = addfile (dirm+film,"r")
115 ;---------------------------------------------------
116 ; Units for these variables are:
126 nsec_per_year = 60*60*24*365
128 ; change unit to g C/m^2/year
130 NPP1 = NPP1 * nsec_per_year
131 NPP2 = NPP2 * nsec_per_year
132 NEE2 = NEE2 * nsec_per_year
136 ;*******************************************************************
137 ; Calculate "nice" bins for binning the data in equally spaced ranges
138 ;********************************************************************
140 ; using model biome class
143 range = fspan(0,nclass,nclass+1)
146 ; Use this range information to grab all the values in a
147 ; particular range, and then take an average.
149 nx = dimsizes(range) - 1
151 ;==============================
153 ;==============================
155 ; using observed biome class
156 ; base = ndtooned(classob)
157 ; using model biome class
158 base = ndtooned(classmod)
160 area_1d = ndtooned(area)
164 yvalues = new((/data_n,nx/),float)
165 yvalues_t = new((/data_n,nx/),float)
168 ; Loop through each range, using base.
172 if (i.ne.(nx-1)) then
173 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
175 idx = ind(base.ge.range(i))
181 data = ndtooned(area)
185 data = ndtooned(NPP1)
189 data = ndtooned(VegC)
193 data = ndtooned(LiCwC)
197 data = ndtooned(SoilC)
201 data = ndtooned(NPP2)
205 data = ndtooned(NEE2)
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 ;#############################################################
252 ;----------------------------------------------------------------
255 good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
258 area_g = yvalues(0,good)
259 NPP1_g = yvalues(1,good)
260 VegC_g = yvalues(2,good)
261 LiCwC_g = yvalues(3,good)
262 SoilC_g = yvalues(4,good)
263 NPP2_g = yvalues(5,good)
264 NEE2_g = yvalues(6,good)
266 n_biome = dimsizes(NPP1_g)
268 NPP_ratio = NPP2_g/NPP1_g
270 ;-----------------------------------------------------------------
273 ; change unit from g to Pg (Peta gram)
276 NPP1_t = yvalues_t(1,good) * factor_unit
277 VegC_t = yvalues_t(2,good) * factor_unit
278 LiCwC_t = yvalues_t(3,good) * factor_unit
279 SoilC_t = yvalues_t(4,good) * factor_unit
280 NEE2_t = yvalues_t(6,good) * factor_unit
285 ;-------------------------------------------------------------
288 ; column (not including header column)
290 col_head = (/"Area (1.e12m2)" \
293 ,"Litter+CWD (gC/m2)" \
299 ncol = dimsizes(col_head)
301 ; row (not including header row)
303 ; using model biome class:
304 row_head = (/"Not Vegetated" \
305 ,"Needleleaf Evergreen Temperate Tree" \
306 ,"Needleleaf Evergreen Boreal Tree" \
307 ; ,"Needleleaf Deciduous Boreal Tree" \
308 ,"Broadleaf Evergreen Tropical Tree" \
309 ,"Broadleaf Evergreen Temperate Tree" \
310 ,"Broadleaf Deciduous Tropical Tree" \
311 ,"Broadleaf Deciduous Temperate Tree" \
312 ; ,"Broadleaf Deciduous Boreal Tree" \
313 ; ,"Broadleaf Evergreen Shrub" \
314 ,"Broadleaf Deciduous Temperate Shrub" \
315 ,"Broadleaf Deciduous Boreal Shrub" \
317 ,"C3 Non-Arctic Grass" \
323 nrow = dimsizes(row_head)
325 ; arrays to be passed to table.
326 text4 = new ((/nrow, ncol/),string )
329 text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)
330 text4(i,1) = sprintf("%.1f",NPP1_g(i))
331 text4(i,2) = sprintf("%.1f",VegC_g(i))
332 text4(i,3) = sprintf("%.1f",LiCwC_g(i))
333 text4(i,4) = sprintf("%.1f",SoilC_g(i))
334 text4(i,5) = sprintf("%.2f",NPP_ratio(i))
335 text4(i,6) = sprintf("%.1f",NEE2_g(i))
338 ;-------------------------------------------------------
341 header_text = "<H1>NEE and Carbon Stocks and Fluxes: Model "+model_name+"</H1>"
343 header = (/"<HTML>" \
345 ,"<TITLE>CLAMP metrics</TITLE>" \
352 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
354 ," <th bgcolor=DDDDDD >Biome Type</th>" \
355 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
356 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
357 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
358 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
359 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
360 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
361 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
364 table_footer = "</table>"
368 lines = new(50000,string)
371 set_line(lines,nline,header)
372 set_line(lines,nline,table_header)
374 ;----------------------------
378 set_line(lines,nline,row_header)
389 set_line(lines,nline,"<th>"+txt0+"</th>")
390 set_line(lines,nline,"<th>"+txt1+"</th>")
391 set_line(lines,nline,"<th>"+txt2+"</th>")
392 set_line(lines,nline,"<th>"+txt3+"</th>")
393 set_line(lines,nline,"<th>"+txt4+"</th>")
394 set_line(lines,nline,"<th>"+txt5+"</th>")
395 set_line(lines,nline,"<th>"+txt6+"</th>")
396 set_line(lines,nline,"<th>"+txt7+"</th>")
398 set_line(lines,nline,row_footer)
400 ;----------------------------
401 set_line(lines,nline,table_footer)
402 set_line(lines,nline,footer)
404 ; Now write to an HTML file.
406 output_html = "table_carbon_sink1.html"
408 idx = ind(.not.ismissing(lines))
409 if(.not.any(ismissing(idx))) then
410 asciiwrite(output_html,lines(idx))
420 delete (table_header)
422 ;-----------------------------------------------------------------
425 ; column (not including header column)
427 col_head = (/"NPP (PgC/yr)" \
429 ,"Litter+CWD (PgC)" \
437 ncol = dimsizes(col_head)
439 ; row (not including header row)
441 ; using model biome class:
442 row_head = (/"Not Vegetated" \
443 ,"Needleleaf Evergreen Temperate Tree" \
444 ,"Needleleaf Evergreen Boreal Tree" \
445 ; ,"Needleleaf Deciduous Boreal Tree" \
446 ,"Broadleaf Evergreen Tropical Tree" \
447 ,"Broadleaf Evergreen Temperate Tree" \
448 ,"Broadleaf Deciduous Tropical Tree" \
449 ,"Broadleaf Deciduous Temperate Tree" \
450 ; ,"Broadleaf Deciduous Boreal Tree" \
451 ; ,"Broadleaf Evergreen Shrub" \
452 ,"Broadleaf Deciduous Temperate Shrub" \
453 ,"Broadleaf Deciduous Boreal Shrub" \
455 ,"C3 Non-Arctic Grass" \
461 nrow = dimsizes(row_head)
463 ; arrays to be passed to table.
464 text4 = new ((/nrow, ncol/),string )
467 text4(i,0) = sprintf("%.1f",NPP1_t(i))
468 text4(i,1) = sprintf("%.1f",VegC_t(i))
469 text4(i,2) = sprintf("%.1f",LiCwC_t(i))
470 text4(i,3) = sprintf("%.1f",SoilC_t(i))
471 text4(i,4) = sprintf("%.1f",NEE2_t(i))
472 text4(i,5) = "<a href=./NPP_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NPP_annual_biome_"+i+".png>annual_plot</a>"
473 text4(i,6) = "<a href=./NEE_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NEE_annual_biome_"+i+".png>annual_plot</a>"
474 text4(i,7) = "<a href=./Fire_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./Fire_annual_biome_"+i+".png>annual_plot</a>"
476 text4(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))
477 text4(nrow-1,1) = sprintf("%.1f",sum(VegC_t))
478 text4(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))
479 text4(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))
480 text4(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))
481 text4(nrow-1,5) = "<a href=./NPP_monthly_global.png>monthly_plot</a> <br> <a href=./NPP_annual_global.png>annual_plot</a>"
482 text4(nrow-1,6) = "<a href=./NEE_monthly_global.png>monthly_plot</a> <br> <a href=./NEE_annual_global.png>annual_plot</a>"
483 text4(nrow-1,7) = "<a href=./Fire_monthly_global.png>monthly_plot</a> <br> <a href=./Fire_annual_global.png>annual_plot</a>"
485 ;**************************************************
487 ;**************************************************
489 header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome): Model "+model_name+"</H1>"
491 header = (/"<HTML>" \
493 ,"<TITLE>CLAMP metrics</TITLE>" \
500 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
502 ," <th bgcolor=DDDDDD >Biome Type</th>" \
503 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
504 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
505 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
506 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
507 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
508 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
509 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
510 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
513 table_footer = "</table>"
517 lines = new(50000,string)
520 set_line(lines,nline,header)
521 set_line(lines,nline,table_header)
522 ;-----------------------------------------------
526 set_line(lines,nline,row_header)
538 set_line(lines,nline,"<th>"+txt0+"</th>")
539 set_line(lines,nline,"<th>"+txt1+"</th>")
540 set_line(lines,nline,"<th>"+txt2+"</th>")
541 set_line(lines,nline,"<th>"+txt3+"</th>")
542 set_line(lines,nline,"<th>"+txt4+"</th>")
543 set_line(lines,nline,"<th>"+txt5+"</th>")
544 set_line(lines,nline,"<th>"+txt6+"</th>")
545 set_line(lines,nline,"<th>"+txt7+"</th>")
546 set_line(lines,nline,"<th>"+txt8+"</th>")
548 set_line(lines,nline,row_footer)
550 ;-----------------------------------------------
551 set_line(lines,nline,table_footer)
552 set_line(lines,nline,footer)
554 ; Now write to an HTML file.
556 output_html = "table_carbon_sink2.html"
558 idx = ind(.not.ismissing(lines))
559 if(.not.any(ismissing(idx))) then
560 asciiwrite(output_html,lines(idx))
567 ;---------------------------------------------------
568 ; read data: model, time series
570 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
571 film = model_name2 + "_Fire_C_1979-2004_monthly.nc"
572 fm = addfile (dirm+film,"r")
576 Fire = fm->COL_FIRE_CLOSS
580 ; Units for these variables are:
586 nsec_per_month = 60*60*24*30
588 ; change unit to g C/m^2/month
590 NPP3 = NPP3 * nsec_per_month
591 NEE3 = NEE3 * nsec_per_month
592 Fire = Fire * nsec_per_month
596 dsizes = dimsizes(NPP3)
599 ntime = nyear * nmonth
604 ;*******************************************************************
605 ; Calculate "nice" bins for binning the data in equally spaced ranges
606 ;********************************************************************
608 ; using model biome class
611 range = fspan(0,nclass,nclass+1)
614 ; Use this range information to grab all the values in a
615 ; particular range, and then take an average.
617 nx = dimsizes(range) - 1
619 ;==============================
621 ;==============================
623 ; using observed biome class
624 ; base = ndtooned(classob)
625 ; using model biome class
626 base = ndtooned(classmod)
630 yvalues = new((/ntime,data_n,nx/),float)
632 ; Loop through each range, using base.
636 if (i.ne.(nx-1)) then
637 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
639 idx = ind(base.ge.range(i))
651 data = ndtooned(NPP3(m,k,:,:))
655 data = ndtooned(NEE3(m,k,:,:))
659 data = ndtooned(Fire(m,k,:,:))
664 if (.not.any(ismissing(idx))) then
665 yvalues(t,n,i) = sum(data(idx)*area_1d(idx))
667 yvalues(t,n,i) = yvalues@_FillValue
670 ;#############################################################
671 ; using model biome class:
672 ; set the following 4 classes to _FillValue:
673 ; (3)Needleleaf Deciduous Boreal Tree,
674 ; (8)Broadleaf Deciduous Boreal Tree,
675 ; (9)Broadleaf Evergreen Shrub,
678 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
679 yvalues(t,n,i) = yvalues@_FillValue
681 ;#############################################################
697 ;----------------------------------------------------------------
698 ; data for tseries plot
700 yvalues_g = new((/ntime,data_n,n_biome/),float)
702 yvalues_g@units = "TgC/month"
704 ; change unit to Tg C/month
705 ; change unit from g to Tg (Tera gram)
708 yvalues_g = yvalues(:,:,good) * factor_unit
710 ;*******************************************************************
711 ; general settings for line plot
712 ;*******************************************************************
716 res@xyDashPatterns = (/0/) ; make lines solid
717 res@xyLineThicknesses = (/2.0/) ; make lines thicker
718 res@xyLineColors = (/"blue"/) ; line color
720 res@trXMinF = year_start
721 res@trXMaxF = year_end + 1
723 res@vpHeightF = 0.4 ; change aspect ratio of plot
727 ; res@gsnMaximize = True
729 ;*******************************************************************
730 ; (A) 1 component in each biome: monthly
731 ;*******************************************************************
733 component = (/"NPP","NEE","Fire"/)
735 ; for x-axis in xyplot
737 timeI = new((/ntime/),integer)
738 timeF = new((/ntime/),float)
739 timeI = ispan(1,ntime,1)
740 timeF = year_start + (timeI-1)/12.
741 timeF@long_name = "year"
743 plot_data = new((/ntime/),float)
744 plot_data@long_name = "TgC/month"
749 plot_name = component(n)+"_monthly_biome_"+ m
751 wks = gsn_open_wks (plot_type,plot_name)
753 title = component(n)+ ": "+ row_head(m)
754 res@tiMainString = title
755 res@tiMainFontHeightF = 0.025
757 plot_data(:) = yvalues_g(:,n,m)
759 plot=gsn_csm_xy(wks,timeF,plot_data,res)
761 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
762 "rm "+plot_name+"."+plot_type)
771 plot_name = component(n)+"_monthly_global"
773 wks = gsn_open_wks (plot_type,plot_name)
775 title = component(n)+ ": Global"
776 res@tiMainString = title
777 res@tiMainFontHeightF = 0.025
780 plot_data(k) = sum(yvalues_g(k,n,:))
783 plot=gsn_csm_xy(wks,timeF,plot_data,res)
785 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
786 "rm "+plot_name+"."+plot_type)
796 ;*******************************************************************
797 ; (B) 1 component in each biome: annually
798 ;*******************************************************************
800 yvalues_a = new((/nyear,data_n,n_biome/),float)
803 yvalues_g!2 = "record"
805 yvalues_a = month_to_annual(yvalues_g,0)
809 ; for x-axis in xyplot
811 timeI = new((/nyear/),integer)
812 timeF = new((/nyear/),float)
813 timeI = ispan(1,nyear,1)
814 timeF = year_start + (timeI-1)
815 timeF@long_name = "year"
817 plot_data = new((/nyear/),float)
818 plot_data@long_name = "TgC/year"
823 plot_name = component(n)+"_annual_biome_"+ m
825 wks = gsn_open_wks (plot_type,plot_name)
827 title = component(n)+ ": "+ row_head(m)
828 res@tiMainString = title
829 res@tiMainFontHeightF = 0.025
831 plot_data(:) = yvalues_a(:,n,m)
833 plot=gsn_csm_xy(wks,timeF,plot_data,res)
835 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
836 "rm "+plot_name+"."+plot_type)
845 plot_name = component(n)+"_annual_global"
847 wks = gsn_open_wks (plot_type,plot_name)
849 title = component(n)+ ": Global"
850 res@tiMainString = title
851 res@tiMainFontHeightF = 0.025
854 plot_data(k) = sum(yvalues_a(k,n,:))
857 plot=gsn_csm_xy(wks,timeF,plot_data,res)
859 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
860 "rm "+plot_name+"."+plot_type)