1 ;********************************************************
3 ; using model biome vlass
5 ; required command line input parameters:
6 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
8 ; histogram normalized by rain and compute correleration
9 ;**************************************************************
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
12 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
13 ;**************************************************************
14 procedure set_line(lines:string,nline:integer,newlines:string)
16 ; add line to ascci/html file
18 nnewlines = dimsizes(newlines)
19 if(nline+nnewlines-1.ge.dimsizes(lines))
20 print("set_line: bad index, not setting anything.")
23 lines(nline:nline+nnewlines-1) = newlines
24 ; print ("lines = " + lines(nline:nline+nnewlines-1))
25 nline = nline + nnewlines
28 ;**************************************************************
35 ;************************************************
37 ;************************************************
42 model_name1 = "i01.06casa"
43 model_name2 = "i01.10casa"
45 ;---------------------------------------------------
46 ; get biome data: model
48 biome_name_mod = "Model PFT Class"
50 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
51 film = "class_pft_"+model_grid+".nc"
53 fm = addfile(dirm+film,"r")
55 classmod = fm->CLASS_PFT
59 ; model data has 17 land-type classes
63 ;--------------------------------------------------
64 ; get model data: landfrac and area
66 dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/"
68 fm_l = addfile (dirm_l+film_l,"r")
70 landfrac = fm_l->landfrac
73 ; change area from km**2 to m**2
76 ;---------------------------------------------------
77 ; take into account landfrac
79 area = area * landfrac
83 ;---------------------------------------------------
84 ; read data: model, group 1
86 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
87 film = model_name1 + "_1980-2004_ANN_climo.nc"
88 fm = addfile (dirm+film,"r")
96 VegC = leafc + woodc + frootc
101 LiCwC = litterc + cwdc
106 ;---------------------------------------------------
107 ; read data: model, group 2
109 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
110 film = model_name2 + "_1990-2004_ANN_climo.nc"
111 fm = addfile (dirm+film,"r")
116 ;---------------------------------------------------
117 ; Units for these variables are:
127 nsec_per_year = 60*60*24*365
129 ; change unit to g C/m^2/year
131 NPP1 = NPP1 * nsec_per_year
132 NPP2 = NPP2 * nsec_per_year
133 NEE2 = NEE2 * nsec_per_year
137 ;*******************************************************************
138 ; Calculate "nice" bins for binning the data in equally spaced ranges
139 ;********************************************************************
141 ; using model biome class
144 range = fspan(0,nclass,nclass+1)
147 ; Use this range information to grab all the values in a
148 ; particular range, and then take an average.
150 nx = dimsizes(range) - 1
152 ;==============================
154 ;==============================
156 ; using observed biome class
157 ; base = ndtooned(classob)
158 ; using model biome class
159 base = ndtooned(classmod)
161 area_1d = ndtooned(area)
165 yvalues = new((/data_n,nx/),float)
166 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>"
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) = "--"
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
597 dsizes = dimsizes(NPP3)
600 ntime = nyear * nmonth
605 ;*******************************************************************
606 ; Calculate "nice" bins for binning the data in equally spaced ranges
607 ;********************************************************************
609 ; using model biome class
612 range = fspan(0,nclass,nclass+1)
615 ; Use this range information to grab all the values in a
616 ; particular range, and then take an average.
618 nx = dimsizes(range) - 1
620 ;==============================
622 ;==============================
624 ; using observed biome class
625 ; base = ndtooned(classob)
626 ; using model biome class
627 base = ndtooned(classmod)
631 yvalues = new((/ntime,data_n,nx/),float)
633 ; Loop through each range, using base.
637 if (i.ne.(nx-1)) then
638 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
640 idx = ind(base.ge.range(i))
652 data = ndtooned(NPP3(m,k,:,:))
656 data = ndtooned(NEE3(m,k,:,:))
660 ; data = ndtooned(Fire(m,k,:,:))
665 if (.not.any(ismissing(idx))) then
666 yvalues(t,n,i) = sum(data(idx)*area_1d(idx))
668 yvalues(t,n,i) = yvalues@_FillValue
671 ;#############################################################
672 ; using model biome class:
673 ; set the following 4 classes to _FillValue:
674 ; (3)Needleleaf Deciduous Boreal Tree,
675 ; (8)Broadleaf Deciduous Boreal Tree,
676 ; (9)Broadleaf Evergreen Shrub,
679 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
680 yvalues(t,n,i) = yvalues@_FillValue
682 ;#############################################################
698 ;----------------------------------------------------------------
699 ; data for tseries plot
701 yvalues_g = new((/ntime,data_n,n_biome/),float)
703 yvalues_g@units = "TgC/month"
705 ; change unit to Tg C/month
706 ; change unit from g to Tg (Tera gram)
709 yvalues_g = yvalues(:,:,good) * factor_unit
711 ;*******************************************************************
712 ; general settings for line plot
713 ;*******************************************************************
717 res@xyDashPatterns = (/0/) ; make lines solid
718 res@xyLineThicknesses = (/2.0/) ; make lines thicker
719 res@xyLineColors = (/"blue"/) ; line color
721 res@trXMinF = year_start
722 res@trXMaxF = year_end + 1
724 res@vpHeightF = 0.4 ; change aspect ratio of plot
728 ; res@gsnMaximize = True
730 ;*******************************************************************
731 ; (A) 1 component in each biome: monthly
732 ;*******************************************************************
734 ; component = (/"NPP","NEE","Fire"/)
735 component = (/"NPP","NEE"/)
737 ; for x-axis in xyplot
739 timeI = new((/ntime/),integer)
740 timeF = new((/ntime/),float)
741 timeI = ispan(1,ntime,1)
742 timeF = year_start + (timeI-1)/12.
743 timeF@long_name = "year"
745 plot_data = new((/ntime/),float)
746 plot_data@long_name = "TgC/month"
751 plot_name = component(n)+"_monthly_biome_"+ m
753 wks = gsn_open_wks (plot_type,plot_name)
755 title = component(n)+ ": "+ row_head(m)
756 res@tiMainString = title
757 res@tiMainFontHeightF = 0.025
759 plot_data(:) = yvalues_g(:,n,m)
761 plot=gsn_csm_xy(wks,timeF,plot_data,res)
763 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
764 "rm "+plot_name+"."+plot_type)
773 plot_name = component(n)+"_monthly_global"
775 wks = gsn_open_wks (plot_type,plot_name)
777 title = component(n)+ ": Global"
778 res@tiMainString = title
779 res@tiMainFontHeightF = 0.025
782 plot_data(k) = sum(yvalues_g(k,n,:))
785 plot=gsn_csm_xy(wks,timeF,plot_data,res)
787 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
788 "rm "+plot_name+"."+plot_type)
798 ;*******************************************************************
799 ; (B) 1 component in each biome: annually
800 ;*******************************************************************
802 yvalues_a = new((/nyear,data_n,n_biome/),float)
805 yvalues_g!2 = "record"
807 yvalues_a = month_to_annual(yvalues_g,0)
811 ; for x-axis in xyplot
813 timeI = new((/nyear/),integer)
814 timeF = new((/nyear/),float)
815 timeI = ispan(1,nyear,1)
816 timeF = year_start + (timeI-1)
817 timeF@long_name = "year"
819 plot_data = new((/nyear/),float)
820 plot_data@long_name = "TgC/year"
825 plot_name = component(n)+"_annual_biome_"+ m
827 wks = gsn_open_wks (plot_type,plot_name)
829 title = component(n)+ ": "+ row_head(m)
830 res@tiMainString = title
831 res@tiMainFontHeightF = 0.025
833 plot_data(:) = yvalues_a(:,n,m)
835 plot=gsn_csm_xy(wks,timeF,plot_data,res)
837 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
838 "rm "+plot_name+"."+plot_type)
847 plot_name = component(n)+"_annual_global"
849 wks = gsn_open_wks (plot_type,plot_name)
851 title = component(n)+ ": Global"
852 res@tiMainString = title
853 res@tiMainFontHeightF = 0.025
856 plot_data(k) = sum(yvalues_a(k,n,:))
859 plot=gsn_csm_xy(wks,timeF,plot_data,res)
861 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
862 "rm "+plot_name+"."+plot_type)