diff -r 000000000000 -r 0c6405ab2ff4 carbon_sink/20x.table+tseries.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/carbon_sink/20x.table+tseries.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,869 @@ +;******************************************************** +; for casa (no fire) +; using model biome vlass +; +; required command line input parameters: +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl +; +; histogram normalized by rain and compute correleration +;************************************************************** +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" + +;************************************************ +; model name and grid +;************************************************ + + model_grid = "T42" + + model_name = "casa" + model_name1 = "i01.06casa" + model_name2 = "i01.10casa" + +;--------------------------------------------------- +; get biome data: model + + biome_name_mod = "Model PFT Class" + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = "class_pft_"+model_grid+".nc" + + fm = addfile(dirm+film,"r") + + classmod = fm->CLASS_PFT + + delete (fm) + +; model data has 17 land-type classes + + nclass_mod = 17 + +;-------------------------------------------------- +; get model data: landfrac and area + + dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/" + film_l = "lnd_T42.nc" + fm_l = addfile (dirm_l+film_l,"r") + + landfrac = fm_l->landfrac + area = fm_l->area + +; 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 + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = model_name1 + "_1980-2004_ANN_climo.nc" + fm = addfile (dirm+film,"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 + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = model_name2 + "_1990-2004_ANN_climo.nc" + fm = addfile (dirm+film,"r") + + NPP2 = fm->NPP + NEE2 = fm->NEE + +;--------------------------------------------------- +; Units for these variables are: + +;NPP1: g C/m^2/s +;NPP2: g C/m^2/s +;NEE2: 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 + NPP2 = NPP2 * nsec_per_year + NEE2 = NEE2 * nsec_per_year + + data_n = 7 + +;******************************************************************* +; 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) + yvalues_t = new((/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 + + 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 + +; 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) + +;---------------------------------------------------------------- +; 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) + + n_biome = dimsizes(NPP1_g) + + NPP_ratio = NPP2_g/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 + + 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)" \ + /) + + 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. + text4 = new ((/nrow, ncol/),string ) + + do i=0,nrow-2 + text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12) + text4(i,1) = sprintf("%.1f",NPP1_g(i)) + text4(i,2) = sprintf("%.1f",VegC_g(i)) + text4(i,3) = sprintf("%.1f",LiCwC_g(i)) + text4(i,4) = sprintf("%.1f",SoilC_g(i)) + text4(i,5) = sprintf("%.2f",NPP_ratio(i)) + text4(i,6) = sprintf("%.1f",NEE2_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)+"
" + 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 = text4(n,0) + txt2 = text4(n,1) + txt3 = text4(n,2) + txt4 = text4(n,3) + txt5 = text4(n,4) + txt6 = text4(n,5) + txt7 = text4(n,6) + + 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,row_footer) + end do +;---------------------------- + set_line(lines,nline,table_footer) + set_line(lines,nline,footer) + +; Now write to an HTML file. + + output_html = "table_carbon_sink1.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 (text4) + 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)" \ + ,"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. + text4 = new ((/nrow, ncol/),string ) + + do i=0,nrow-2 + text4(i,0) = sprintf("%.1f",NPP1_t(i)) + text4(i,1) = sprintf("%.1f",VegC_t(i)) + text4(i,2) = sprintf("%.1f",LiCwC_t(i)) + text4(i,3) = sprintf("%.1f",SoilC_t(i)) + text4(i,4) = sprintf("%.1f",NEE2_t(i)) + text4(i,5) = "monthly_plot
annual_plot" + text4(i,6) = "monthly_plot
annual_plot" + text4(i,7) = "--" + end do + text4(nrow-1,0) = sprintf("%.1f",sum(NPP1_t)) + text4(nrow-1,1) = sprintf("%.1f",sum(VegC_t)) + text4(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t)) + text4(nrow-1,3) = sprintf("%.1f",sum(SoilC_t)) + text4(nrow-1,4) = sprintf("%.1f",sum(NEE2_t)) + text4(nrow-1,5) = "monthly_plot
annual_plot" + text4(nrow-1,6) = "monthly_plot
annual_plot" + text4(nrow-1,7) = "--" + +;************************************************** +; 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)+"
" + 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 = text4(n,0) + txt2 = text4(n,1) + txt3 = text4(n,2) + txt4 = text4(n,3) + txt5 = text4(n,4) + txt6 = text4(n,5) + txt7 = text4(n,6) + txt8 = text4(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_carbon_sink2.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 data: model, time series + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = model_name2 + "_Fire_C_1979-2004_monthly.nc" + fm = addfile (dirm+film,"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 + NEE3 = NEE3 * nsec_per_month +;Fire = Fire * nsec_per_month + +;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) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + + clear (wks) + delete (plot) + 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) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + + clear (wks) + delete (plot) + 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) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + + clear (wks) + delete (plot) + 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) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ + "rm "+plot_name+"."+plot_type) + + clear (wks) + delete (plot) + end do + +end +