diff -r 000000000000 -r 0c6405ab2ff4 all/09x.carbon_sink.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/all/09x.carbon_sink.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,890 @@ +;************************************************************** +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" +;************************************************************** +procedure set_line(lines:string,nline:integer,newlines:string) +begin +; add line to ascci/html file + + nnewlines = dimsizes(newlines) + if(nline+nnewlines-1.ge.dimsizes(lines)) + print("set_line: bad index, not setting anything.") + return + end if + lines(nline:nline+nnewlines-1) = newlines +; print ("lines = " + lines(nline:nline+nnewlines-1)) + nline = nline + nnewlines + return +end +;************************************************************** +; Main code. +begin + + plot_type = "ps" + plot_type_new = "png" + +;---------------------------------------------------------- +; edit current model for movel1_vs_model2 + + if (isvar("compare")) then + html_name2 = compare+"/table.html" + html_new2 = html_name2 +".new" + end if + +;---------------------------------------------------------- +; edit table.html for current model + + html_name = model_name+"/table.html" + html_new = html_name +".new" + +;---------------------------------------------------------- +; get biome data: model + + biome_name_mod = "Model PFT Class" + + film_c = "class_pft_"+ model_grid +".nc" + fm_c = addfile (dirs+film_c,"r") + classmod = fm_c->CLASS_PFT + + delete (fm_c) + +; model data has 17 land-type classes + nclass_mod = 17 + +;------------------------------------ +; get model landfrac and area + + film_l = "lnd_"+ model_grid +".nc" + fm_l = addfile (dirs+film_l,"r") + landfrac = fm_l->landfrac + area = fm_l->area + + delete (fm_l) + +; change area from km**2 to m**2 + area = area * 1.e6 + +;------------------------------------- +; take into account landfrac + +; area = area * landfrac + +; delete (landfrac) + +;--------------------------------------------------- +; read data: model, group 1 + + fm = addfile (dirm+film4,"r") + + NPP1 = fm->NPP + + leafc = fm->LEAFC + woodc = fm->WOODC + frootc = fm->FROOTC + VegC = leafc + VegC = leafc + woodc + frootc + + litterc = fm->LITTERC + cwdc = fm->CWDC + LiCwC = litterc + LiCwC = litterc + cwdc + + SoilC = fm->SOILC + + delete (fm) +;--------------------------------------------------- +; read data: model, group 2 + + fm = addfile (dirm+film5,"r") + + NPP2 = fm->NPP + NEE2 = fm->NEE + GPP2 = fm->GPP + + delete (fm) +;--------------------------------------------------- +; Units for these variables are: + +;NPP1: g C/m^2/s +;NPP2: g C/m^2/s +;NEE2: g C/m^2/s +;GPP2: g C/m^2/s + +;VegC: g C/m^2 +;LiCwC: g C/m^2 +;SoilC: g C/m^2 + + nsec_per_year = 60*60*24*365 + +; change unit to g C/m^2/year + + NPP1 = NPP1 * nsec_per_year * conform(NPP1,landfrac,(/1,2/)) + NPP2 = NPP2 * nsec_per_year * conform(NPP2,landfrac,(/1,2/)) + NEE2 = NEE2 * nsec_per_year * conform(NEE2,landfrac,(/1,2/)) + GPP2 = GPP2 * nsec_per_year * conform(GPP2,landfrac,(/1,2/)) + + VegC = VegC * conform(VegC,landfrac,(/1,2/)) + LiCwC = LiCwC * conform(LiCwC,landfrac,(/1,2/)) + SoilC = SoilC * conform(SoilC,landfrac,(/1,2/)) + + data_n = 8 + +;******************************************************************* +; Calculate "nice" bins for binning the data in equally spaced ranges +;******************************************************************** + +; using model biome class + nclass = nclass_mod + + range = fspan(0,nclass,nclass+1) + +; print (range) +; Use this range information to grab all the values in a +; particular range, and then take an average. + + nx = dimsizes(range) - 1 + +;============================== +; put data into bins +;============================== + +; using observed biome class +; base = ndtooned(classob) +; using model biome class + base = ndtooned(classmod) + + area_1d = ndtooned(area) + +; output + + yvalues = new((/data_n,nx/),float) ; (per m2) + yvalues_t = new((/data_n,nx/),float) ; (per biome) + +; Loop through each range, using base. + + do i=0,nx-1 + + if (i.ne.(nx-1)) then + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) + else + idx = ind(base.ge.range(i)) + end if + + do n = 0,data_n-1 + + if (n.eq.0) then + data = ndtooned(area) + end if + + if (n.eq.1) then + data = ndtooned(NPP1) + end if + + if (n.eq.2) then + data = ndtooned(VegC) + end if + + if (n.eq.3) then + data = ndtooned(LiCwC) + end if + + if (n.eq.4) then + data = ndtooned(SoilC) + end if + + if (n.eq.5) then + data = ndtooned(NPP2) + end if + + if (n.eq.6) then + data = ndtooned(NEE2) + end if + + if (n.eq.7) then + data = ndtooned(GPP2) + end if + +; Calculate sum and average + + if (.not.any(ismissing(idx))) then + if (n.eq.0) then + yvalues(n,i) = sum(data(idx)) + yvalues_t(n,i) = sum(data(idx)) + else + yvalues(n,i) = avg(data(idx)) + yvalues_t(n,i) = sum(data(idx)*area_1d(idx)) + end if + else + yvalues(n,i) = yvalues@_FillValue + yvalues_t(n,i) = yvalues@_FillValue + end if + +;############################################################# +; using model biome class: +; set the following 4 classes to _FillValue: +; (3)Needleleaf Deciduous Boreal Tree, +; (8)Broadleaf Deciduous Boreal Tree, +; (9)Broadleaf Evergreen Shrub, +; (16)Wheat + + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then + yvalues(n,i) = yvalues@_FillValue + yvalues_t(n,i) = yvalues@_FillValue + end if +;############################################################# + + delete (data) + end do + + delete (idx) + end do + + delete (base) + delete (area) + delete (NPP1) + delete (VegC) + delete (LiCwC) + delete (SoilC) + delete (NPP2) + delete (NEE2) + delete (GPP2) + +;---------------------------------------------------------------- +; data for table1 + + good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:))) +;print (good) + + area_g = yvalues(0,good) + NPP1_g = yvalues(1,good) + VegC_g = yvalues(2,good) + LiCwC_g = yvalues(3,good) + SoilC_g = yvalues(4,good) + NPP2_g = yvalues(5,good) + NEE2_g = yvalues(6,good) + GPP2_g = yvalues(7,good) + + NPP_ratio = NPP2_g/NPP1_g + + n_biome = dimsizes(NPP1_g) + +;----------------------------------------------------------------- +; data for table2 + +; change unit from g to Pg (Peta gram) + factor_unit = 1.e-15 + + NPP1_t = yvalues_t(1,good) * factor_unit + VegC_t = yvalues_t(2,good) * factor_unit + LiCwC_t = yvalues_t(3,good) * factor_unit + SoilC_t = yvalues_t(4,good) * factor_unit + NEE2_t = yvalues_t(6,good) * factor_unit + GPP2_t = yvalues_t(7,good) * factor_unit + + delete (yvalues) + delete (yvalues_t) + +;------------------------------------------------------------- +; html table1 data + +; column (not including header column) + + col_head = (/"Area (1.e12m2)" \ + ,"NPP (gC/m2/yr)" \ + ,"VegC (gC/m2)" \ + ,"Litter+CWD (gC/m2)" \ + ,"SoilC (gC/m2)" \ + ,"NPP_ratio" \ + ,"NEE (gC/m2/yr)" \ + ,"GPP (gC/m2/yr)" \ + /) + + ncol = dimsizes(col_head) + +; row (not including header row) + +; using model biome class: + row_head = (/"Not Vegetated" \ + ,"Needleleaf Evergreen Temperate Tree" \ + ,"Needleleaf Evergreen Boreal Tree" \ +; ,"Needleleaf Deciduous Boreal Tree" \ + ,"Broadleaf Evergreen Tropical Tree" \ + ,"Broadleaf Evergreen Temperate Tree" \ + ,"Broadleaf Deciduous Tropical Tree" \ + ,"Broadleaf Deciduous Temperate Tree" \ +; ,"Broadleaf Deciduous Boreal Tree" \ +; ,"Broadleaf Evergreen Shrub" \ + ,"Broadleaf Deciduous Temperate Shrub" \ + ,"Broadleaf Deciduous Boreal Shrub" \ + ,"C3 Arctic Grass" \ + ,"C3 Non-Arctic Grass" \ + ,"C4 Grass" \ + ,"Corn" \ +; ,"Wheat" \ + ,"All Biome" \ + /) + nrow = dimsizes(row_head) + +; arrays to be passed to table. + text = new ((/nrow, ncol/),string ) + + do i=0,nrow-2 + text(i,0) = sprintf("%.1f",area_g(i)*1.e-12) + text(i,1) = sprintf("%.1f",NPP1_g(i)) + text(i,2) = sprintf("%.1f",VegC_g(i)) + text(i,3) = sprintf("%.1f",LiCwC_g(i)) + text(i,4) = sprintf("%.1f",SoilC_g(i)) + text(i,5) = sprintf("%.2f",NPP_ratio(i)) + text(i,6) = sprintf("%.1f",NEE2_g(i)) + text(i,7) = sprintf("%.1f",GPP2_g(i)) + end do + +;------------------------------------------------------- +; create html table1 + + header_text = "

NEE and Carbon Stocks and Fluxes: Model "+model_name+"

" + + header = (/"" \ + ,"" \ + ,"CLAMP metrics" \ + ,"" \ + ,header_text \ + /) + footer = "" + + table_header = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
Biome Type"+col_head(0)+""+col_head(1)+""+col_head(2)+""+col_head(3)+""+col_head(4)+""+col_head(5)+""+col_head(6)+""+col_head(7)+"
" + row_header = "" + row_footer = "" + + lines = new(50000,string) + nline = 0 + + set_line(lines,nline,header) + set_line(lines,nline,table_header) + +;---------------------------- +;row of table + + do n = 0,nrow-2 + set_line(lines,nline,row_header) + + txt0 = row_head(n) + txt1 = text(n,0) + txt2 = text(n,1) + txt3 = text(n,2) + txt4 = text(n,3) + txt5 = text(n,4) + txt6 = text(n,5) + txt7 = text(n,6) + txt8 = text(n,7) + + set_line(lines,nline,""+txt0+"") + set_line(lines,nline,""+txt1+"") + set_line(lines,nline,""+txt2+"") + set_line(lines,nline,""+txt3+"") + set_line(lines,nline,""+txt4+"") + set_line(lines,nline,""+txt5+"") + set_line(lines,nline,""+txt6+"") + set_line(lines,nline,""+txt7+"") + set_line(lines,nline,""+txt8+"") + + set_line(lines,nline,row_footer) + end do +;---------------------------- + set_line(lines,nline,table_footer) + set_line(lines,nline,footer) + +; Now write to an HTML file. + + output_html = "table_per_m2.html" + + idx = ind(.not.ismissing(lines)) + if(.not.any(ismissing(idx))) then + asciiwrite(output_html,lines(idx)) + else + print ("error?") + end if + + delete (idx) + + delete (col_head) + delete (row_head) + delete (text) + delete (table_header) + +;----------------------------------------------------------------- +; html table2 data + +; column (not including header column) + + col_head = (/"NPP (PgC/yr)" \ + ,"VegC (PgC)" \ + ,"Litter+CWD (PgC)" \ + ,"SoilC (PgC)" \ + ,"NEE (PgC/yr)" \ + ,"GPP (PgC/yr)" \ + ,"NPP timeseries" \ + ,"NEE timeseries" \ + ,"Fire timeseries" \ + /) + + ncol = dimsizes(col_head) + +; row (not including header row) + +; using model biome class: + row_head = (/"Not Vegetated" \ + ,"Needleleaf Evergreen Temperate Tree" \ + ,"Needleleaf Evergreen Boreal Tree" \ +; ,"Needleleaf Deciduous Boreal Tree" \ + ,"Broadleaf Evergreen Tropical Tree" \ + ,"Broadleaf Evergreen Temperate Tree" \ + ,"Broadleaf Deciduous Tropical Tree" \ + ,"Broadleaf Deciduous Temperate Tree" \ +; ,"Broadleaf Deciduous Boreal Tree" \ +; ,"Broadleaf Evergreen Shrub" \ + ,"Broadleaf Deciduous Temperate Shrub" \ + ,"Broadleaf Deciduous Boreal Shrub" \ + ,"C3 Arctic Grass" \ + ,"C3 Non-Arctic Grass" \ + ,"C4 Grass" \ + ,"Corn" \ +; ,"Wheat" \ + ,"All Biome" \ + /) + nrow = dimsizes(row_head) + +; arrays to be passed to table. + text = new ((/nrow, ncol/),string ) + + do i=0,nrow-2 + text(i,0) = sprintf("%.1f",NPP1_t(i)) + text(i,1) = sprintf("%.1f",VegC_t(i)) + text(i,2) = sprintf("%.1f",LiCwC_t(i)) + text(i,3) = sprintf("%.1f",SoilC_t(i)) + text(i,4) = sprintf("%.1f",NEE2_t(i)) + text(i,5) = sprintf("%.1f",GPP2_t(i)) + text(i,6) = "monthly_plot
annual_plot" + text(i,7) = "monthly_plot
annual_plot" + text(i,8) = "--" + end do + + text(nrow-1,0) = sprintf("%.1f",sum(NPP1_t)) + text(nrow-1,1) = sprintf("%.1f",sum(VegC_t)) + text(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t)) + text(nrow-1,3) = sprintf("%.1f",sum(SoilC_t)) + text(nrow-1,4) = sprintf("%.1f",sum(NEE2_t)) + text(nrow-1,5) = sprintf("%.1f",sum(GPP2_t)) + text(nrow-1,6) = "monthly_plot
annual_plot" + text(nrow-1,7) = "monthly_plot
annual_plot" + text(nrow-1,8) = "--" + +;************************************************** +; create html table2 +;************************************************** + + header_text = "

NEE and Carbon Stocks and Fluxes (per biome): Model "+model_name+"

" + + header = (/"" \ + ,"" \ + ,"CLAMP metrics" \ + ,"" \ + ,header_text \ + /) + footer = "" + + table_header = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
Biome Type"+col_head(0)+""+col_head(1)+""+col_head(2)+""+col_head(3)+""+col_head(4)+""+col_head(5)+""+col_head(6)+""+col_head(7)+""+col_head(8)+"
" + row_header = "" + row_footer = "" + + lines = new(50000,string) + nline = 0 + + set_line(lines,nline,header) + set_line(lines,nline,table_header) +;----------------------------------------------- +;row of table + + do n = 0,nrow-1 + set_line(lines,nline,row_header) + + txt0 = row_head(n) + txt1 = text(n,0) + txt2 = text(n,1) + txt3 = text(n,2) + txt4 = text(n,3) + txt5 = text(n,4) + txt6 = text(n,5) + txt7 = text(n,6) + txt8 = text(n,7) + txt9 = text(n,8) + + set_line(lines,nline,""+txt0+"") + set_line(lines,nline,""+txt1+"") + set_line(lines,nline,""+txt2+"") + set_line(lines,nline,""+txt3+"") + set_line(lines,nline,""+txt4+"") + set_line(lines,nline,""+txt5+"") + set_line(lines,nline,""+txt6+"") + set_line(lines,nline,""+txt7+"") + set_line(lines,nline,""+txt8+"") + set_line(lines,nline,""+txt9+"") + + set_line(lines,nline,row_footer) + end do +;----------------------------------------------- + set_line(lines,nline,table_footer) + set_line(lines,nline,footer) + +; Now write to an HTML file. + + output_html = "table_per_biome.html" + + idx = ind(.not.ismissing(lines)) + if(.not.any(ismissing(idx))) then + asciiwrite(output_html,lines(idx)) + else + print ("error?") + end if + + delete (idx) + +;--------------------------------------------------- +; read model data, time series: + + fm = addfile (dirm+film7,"r") + + NPP3 = fm->NPP + NEE3 = fm->NEE +;Fire = fm->COL_FIRE_CLOSS + + delete (fm) + +; Units for these variables are: + +;NPP3: g C/m^2/s +;NEE3: g C/m^2/s +;Fire: g C/m^2/s + + nsec_per_month = 60*60*24*30 + +; change unit to g C/m^2/month + + NPP3 = NPP3 * nsec_per_month * conform(NPP3,landfrac,(/2,3/)) + NEE3 = NEE3 * nsec_per_month * conform(NEE3,landfrac,(/2,3/)) +;Fire = Fire * nsec_per_month * conform(Fire,landfrac,(/2,3/)) + +;data_n = 3 + data_n = 2 + + dsizes = dimsizes(NPP3) + nyear = dsizes(0) + nmonth = dsizes(1) + ntime = nyear * nmonth + + year_start = 1979 + year_end = 2004 + +;******************************************************************* +; Calculate "nice" bins for binning the data in equally spaced ranges +;******************************************************************** + +; using model biome class + nclass = nclass_mod + + range = fspan(0,nclass,nclass+1) + +; print (range) +; Use this range information to grab all the values in a +; particular range, and then take an average. + + nx = dimsizes(range) - 1 + +;============================== +; put data into bins +;============================== + +; using observed biome class +; base = ndtooned(classob) +; using model biome class + base = ndtooned(classmod) + +; output + + yvalues = new((/ntime,data_n,nx/),float) + +; Loop through each range, using base. + + do i=0,nx-1 + + if (i.ne.(nx-1)) then + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) + else + idx = ind(base.ge.range(i)) + end if + + do n = 0,data_n-1 + + t = -1 + do m = 0,nyear-1 + do k = 0,nmonth-1 + + t = t + 1 + + if (n.eq.0) then + data = ndtooned(NPP3(m,k,:,:)) + end if + + if (n.eq.1) then + data = ndtooned(NEE3(m,k,:,:)) + end if + +; if (n.eq.2) then +; data = ndtooned(Fire(m,k,:,:)) +; end if + +; Calculate average + + if (.not.any(ismissing(idx))) then + yvalues(t,n,i) = sum(data(idx)*area_1d(idx)) + else + yvalues(t,n,i) = yvalues@_FillValue + end if + +;############################################################# +; using model biome class: +; set the following 4 classes to _FillValue: +; (3)Needleleaf Deciduous Boreal Tree, +; (8)Broadleaf Deciduous Boreal Tree, +; (9)Broadleaf Evergreen Shrub, +; (16)Wheat + + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then + yvalues(t,n,i) = yvalues@_FillValue + end if +;############################################################# + + end do + end do + + delete(data) + end do + + delete(idx) + end do + + delete (base) + delete (NPP3) + delete (NEE3) +; delete (Fire) + +;---------------------------------------------------------------- +; data for tseries plot + + yvalues_g = new((/ntime,data_n,n_biome/),float) + + yvalues_g@units = "TgC/month" + +; change unit to Tg C/month +; change unit from g to Tg (Tera gram) + factor_unit = 1.e-12 + + yvalues_g = yvalues(:,:,good) * factor_unit + +;******************************************************************* +; general settings for line plot +;******************************************************************* + +; res + res = True + res@xyDashPatterns = (/0/) ; make lines solid + res@xyLineThicknesses = (/2.0/) ; make lines thicker + res@xyLineColors = (/"blue"/) ; line color + + res@trXMinF = year_start + res@trXMaxF = year_end + 1 + + res@vpHeightF = 0.4 ; change aspect ratio of plot +; res@vpWidthF = 0.8 + res@vpWidthF = 0.75 + +; res@gsnMaximize = True + +;******************************************************************* +; (A) 1 component in each biome: monthly +;******************************************************************* + +; component = (/"NPP","NEE","Fire"/) + component = (/"NPP","NEE"/) + +; for x-axis in xyplot + + timeI = new((/ntime/),integer) + timeF = new((/ntime/),float) + timeI = ispan(1,ntime,1) + timeF = year_start + (timeI-1)/12. + timeF@long_name = "year" + + plot_data = new((/ntime/),float) + plot_data@long_name = "TgC/month" + + do n = 0, data_n-1 + do m = 0, n_biome-1 + + plot_name = component(n)+"_monthly_biome_"+ m + + wks = gsn_open_wks (plot_type,plot_name) + + title = component(n)+ ": "+ row_head(m) + res@tiMainString = title + res@tiMainFontHeightF = 0.025 + + plot_data(:) = yvalues_g(:,n,m) + + plot=gsn_csm_xy(wks,timeF,plot_data,res) + + delete (wks) + delete (plot) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + end do + end do + + do n = 0, data_n-1 + + plot_name = component(n)+"_monthly_global" + + wks = gsn_open_wks (plot_type,plot_name) + + title = component(n)+ ": Global" + res@tiMainString = title + res@tiMainFontHeightF = 0.025 + + do k = 0,ntime-1 + plot_data(k) = sum(yvalues_g(k,n,:)) + end do + + plot=gsn_csm_xy(wks,timeF,plot_data,res) + + delete (wks) + delete (plot) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + end do + + delete (plot_data) + delete (timeI) + delete (timeF) + +;******************************************************************* +; (B) 1 component in each biome: annually +;******************************************************************* + + yvalues_a = new((/nyear,data_n,n_biome/),float) + yvalues_g!0 = "time" + yvalues_g!1 = "case" + yvalues_g!2 = "record" + + yvalues_a = month_to_annual(yvalues_g,0) + + delete (yvalues_g) + +; for x-axis in xyplot + + timeI = new((/nyear/),integer) + timeF = new((/nyear/),float) + timeI = ispan(1,nyear,1) + timeF = year_start + (timeI-1) + timeF@long_name = "year" + + plot_data = new((/nyear/),float) + plot_data@long_name = "TgC/year" + + do n = 0, data_n-1 + do m = 0, n_biome-1 + + plot_name = component(n)+"_annual_biome_"+ m + + wks = gsn_open_wks (plot_type,plot_name) + + title = component(n)+ ": "+ row_head(m) + res@tiMainString = title + res@tiMainFontHeightF = 0.025 + + plot_data(:) = yvalues_a(:,n,m) + + plot=gsn_csm_xy(wks,timeF,plot_data,res) + + delete (wks) + delete (plot) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + end do + end do + + do n = 0, data_n-1 + + plot_name = component(n)+"_annual_global" + + wks = gsn_open_wks (plot_type,plot_name) + + title = component(n)+ ": Global" + res@tiMainString = title + res@tiMainFontHeightF = 0.025 + + do k = 0,nyear-1 + plot_data(k) = sum(yvalues_a(k,n,:)) + end do + + plot=gsn_csm_xy(wks,timeF,plot_data,res) + + delete (wks) + delete (plot) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + end do + +;**************************************** +; output plot and html +;**************************************** + output_dir = model_name+"/carbon_sink" + + system("mv *.png *.html " + output_dir) +;**************************************** + +end +