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

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

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

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

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