1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/carbon_sink/20.table+tseries.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,867 @@
1.4 +;********************************************************
1.5 +;using model biome vlass
1.6 +;
1.7 +; required command line input parameters:
1.8 +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
1.9 +;
1.10 +; histogram normalized by rain and compute correleration
1.11 +;**************************************************************
1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.13 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.14 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.15 +;**************************************************************
1.16 +procedure set_line(lines:string,nline:integer,newlines:string)
1.17 +begin
1.18 +; add line to ascci/html file
1.19 +
1.20 + nnewlines = dimsizes(newlines)
1.21 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.22 + print("set_line: bad index, not setting anything.")
1.23 + return
1.24 + end if
1.25 + lines(nline:nline+nnewlines-1) = newlines
1.26 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.27 + nline = nline + nnewlines
1.28 + return
1.29 +end
1.30 +;**************************************************************
1.31 +; Main code.
1.32 +begin
1.33 +
1.34 + plot_type = "ps"
1.35 + plot_type_new = "png"
1.36 +
1.37 +;************************************************
1.38 +; model name and grid
1.39 +;************************************************
1.40 +
1.41 + model_grid = "T42"
1.42 +
1.43 + model_name = "cn"
1.44 + model_name1 = "i01.06cn"
1.45 + model_name2 = "i01.10cn"
1.46 +
1.47 +;---------------------------------------------------
1.48 +; get biome data: model
1.49 +
1.50 + biome_name_mod = "Model PFT Class"
1.51 +
1.52 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.53 + film = "class_pft_"+model_grid+".nc"
1.54 +
1.55 + fm = addfile(dirm+film,"r")
1.56 +
1.57 + classmod = fm->CLASS_PFT
1.58 +
1.59 + delete (fm)
1.60 +
1.61 +; model data has 17 land-type classes
1.62 +
1.63 + nclass_mod = 17
1.64 +
1.65 +;--------------------------------------------------
1.66 +; get model data: landfrac and area
1.67 +
1.68 + dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/"
1.69 + film_l = "lnd_T42.nc"
1.70 + fm_l = addfile (dirm_l+film_l,"r")
1.71 +
1.72 + landfrac = fm_l->landfrac
1.73 + area = fm_l->area
1.74 +
1.75 +; change area from km**2 to m**2
1.76 + area = area * 1.e6
1.77 +
1.78 +;---------------------------------------------------
1.79 +; take into account landfrac
1.80 +
1.81 + area = area * landfrac
1.82 +
1.83 + delete (landfrac)
1.84 +
1.85 +;---------------------------------------------------
1.86 +; read data: model, group 1
1.87 +
1.88 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.89 + film = model_name1 + "_1980-2004_ANN_climo.nc"
1.90 + fm = addfile (dirm+film,"r")
1.91 +
1.92 + NPP1 = fm->NPP
1.93 +
1.94 + leafc = fm->LEAFC
1.95 + woodc = fm->WOODC
1.96 + frootc = fm->FROOTC
1.97 + VegC = leafc
1.98 + VegC = leafc + woodc + frootc
1.99 +
1.100 + litterc = fm->LITTERC
1.101 + cwdc = fm->CWDC
1.102 + LiCwC = litterc
1.103 + LiCwC = litterc + cwdc
1.104 +
1.105 + SoilC = fm->SOILC
1.106 +
1.107 + delete (fm)
1.108 +;---------------------------------------------------
1.109 +; read data: model, group 2
1.110 +
1.111 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.112 + film = model_name2 + "_1990-2004_ANN_climo.nc"
1.113 + fm = addfile (dirm+film,"r")
1.114 +
1.115 + NPP2 = fm->NPP
1.116 + NEE2 = fm->NEE
1.117 +
1.118 +;---------------------------------------------------
1.119 +; Units for these variables are:
1.120 +
1.121 +;NPP1: g C/m^2/s
1.122 +;NPP2: g C/m^2/s
1.123 +;NEE2: g C/m^2/s
1.124 +
1.125 +;VegC: g C/m^2
1.126 +;LiCwC: g C/m^2
1.127 +;SoilC: g C/m^2
1.128 +
1.129 + nsec_per_year = 60*60*24*365
1.130 +
1.131 +; change unit to g C/m^2/year
1.132 +
1.133 + NPP1 = NPP1 * nsec_per_year
1.134 + NPP2 = NPP2 * nsec_per_year
1.135 + NEE2 = NEE2 * nsec_per_year
1.136 +
1.137 + data_n = 7
1.138 +
1.139 +;*******************************************************************
1.140 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.141 +;********************************************************************
1.142 +
1.143 +; using model biome class
1.144 + nclass = nclass_mod
1.145 +
1.146 + range = fspan(0,nclass,nclass+1)
1.147 +
1.148 +; print (range)
1.149 +; Use this range information to grab all the values in a
1.150 +; particular range, and then take an average.
1.151 +
1.152 + nx = dimsizes(range) - 1
1.153 +
1.154 +;==============================
1.155 +; put data into bins
1.156 +;==============================
1.157 +
1.158 +; using observed biome class
1.159 +; base = ndtooned(classob)
1.160 +; using model biome class
1.161 + base = ndtooned(classmod)
1.162 +
1.163 + area_1d = ndtooned(area)
1.164 +
1.165 +; output
1.166 +
1.167 + yvalues = new((/data_n,nx/),float)
1.168 + yvalues_t = new((/data_n,nx/),float)
1.169 +
1.170 +
1.171 +; Loop through each range, using base.
1.172 +
1.173 + do i=0,nx-1
1.174 +
1.175 + if (i.ne.(nx-1)) then
1.176 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.177 + else
1.178 + idx = ind(base.ge.range(i))
1.179 + end if
1.180 +
1.181 + do n = 0,data_n-1
1.182 +
1.183 + if (n.eq.0) then
1.184 + data = ndtooned(area)
1.185 + end if
1.186 +
1.187 + if (n.eq.1) then
1.188 + data = ndtooned(NPP1)
1.189 + end if
1.190 +
1.191 + if (n.eq.2) then
1.192 + data = ndtooned(VegC)
1.193 + end if
1.194 +
1.195 + if (n.eq.3) then
1.196 + data = ndtooned(LiCwC)
1.197 + end if
1.198 +
1.199 + if (n.eq.4) then
1.200 + data = ndtooned(SoilC)
1.201 + end if
1.202 +
1.203 + if (n.eq.5) then
1.204 + data = ndtooned(NPP2)
1.205 + end if
1.206 +
1.207 + if (n.eq.6) then
1.208 + data = ndtooned(NEE2)
1.209 + end if
1.210 +
1.211 +; Calculate sum and average
1.212 +
1.213 + if (.not.any(ismissing(idx))) then
1.214 + if (n.eq.0) then
1.215 + yvalues(n,i) = sum(data(idx))
1.216 + yvalues_t(n,i) = sum(data(idx))
1.217 + else
1.218 + yvalues(n,i) = avg(data(idx))
1.219 + yvalues_t(n,i) = sum(data(idx)*area_1d(idx))
1.220 + end if
1.221 + else
1.222 + yvalues(n,i) = yvalues@_FillValue
1.223 + yvalues_t(n,i) = yvalues@_FillValue
1.224 + end if
1.225 +
1.226 +;#############################################################
1.227 +; using model biome class:
1.228 +; set the following 4 classes to _FillValue:
1.229 +; (3)Needleleaf Deciduous Boreal Tree,
1.230 +; (8)Broadleaf Deciduous Boreal Tree,
1.231 +; (9)Broadleaf Evergreen Shrub,
1.232 +; (16)Wheat
1.233 +
1.234 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.235 + yvalues(n,i) = yvalues@_FillValue
1.236 + yvalues_t(n,i) = yvalues@_FillValue
1.237 + end if
1.238 +;#############################################################
1.239 +
1.240 + delete (data)
1.241 + end do
1.242 +
1.243 + delete (idx)
1.244 + end do
1.245 +
1.246 + delete (base)
1.247 + delete (area)
1.248 + delete (NPP1)
1.249 + delete (VegC)
1.250 + delete (LiCwC)
1.251 + delete (SoilC)
1.252 + delete (NPP2)
1.253 + delete (NEE2)
1.254 +
1.255 +;----------------------------------------------------------------
1.256 +; data for table1
1.257 +
1.258 + good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
1.259 +;print (good)
1.260 +
1.261 + area_g = yvalues(0,good)
1.262 + NPP1_g = yvalues(1,good)
1.263 + VegC_g = yvalues(2,good)
1.264 + LiCwC_g = yvalues(3,good)
1.265 + SoilC_g = yvalues(4,good)
1.266 + NPP2_g = yvalues(5,good)
1.267 + NEE2_g = yvalues(6,good)
1.268 +
1.269 + n_biome = dimsizes(NPP1_g)
1.270 +
1.271 + NPP_ratio = NPP2_g/NPP1_g
1.272 +
1.273 +;-----------------------------------------------------------------
1.274 +; data for table2
1.275 +
1.276 +; change unit from g to Pg (Peta gram)
1.277 + factor_unit = 1.e-15
1.278 +
1.279 + NPP1_t = yvalues_t(1,good) * factor_unit
1.280 + VegC_t = yvalues_t(2,good) * factor_unit
1.281 + LiCwC_t = yvalues_t(3,good) * factor_unit
1.282 + SoilC_t = yvalues_t(4,good) * factor_unit
1.283 + NEE2_t = yvalues_t(6,good) * factor_unit
1.284 +
1.285 + delete (yvalues)
1.286 + delete (yvalues_t)
1.287 +
1.288 +;-------------------------------------------------------------
1.289 +; html table1 data
1.290 +
1.291 +; column (not including header column)
1.292 +
1.293 + col_head = (/"Area (1.e12m2)" \
1.294 + ,"NPP (gC/m2/yr)" \
1.295 + ,"VegC (gC/m2)" \
1.296 + ,"Litter+CWD (gC/m2)" \
1.297 + ,"SoilC (gC/m2)" \
1.298 + ,"NPP_ratio" \
1.299 + ,"NEE (gC/m2/yr)" \
1.300 + /)
1.301 +
1.302 + ncol = dimsizes(col_head)
1.303 +
1.304 +; row (not including header row)
1.305 +
1.306 +; using model biome class:
1.307 + row_head = (/"Not Vegetated" \
1.308 + ,"Needleleaf Evergreen Temperate Tree" \
1.309 + ,"Needleleaf Evergreen Boreal Tree" \
1.310 +; ,"Needleleaf Deciduous Boreal Tree" \
1.311 + ,"Broadleaf Evergreen Tropical Tree" \
1.312 + ,"Broadleaf Evergreen Temperate Tree" \
1.313 + ,"Broadleaf Deciduous Tropical Tree" \
1.314 + ,"Broadleaf Deciduous Temperate Tree" \
1.315 +; ,"Broadleaf Deciduous Boreal Tree" \
1.316 +; ,"Broadleaf Evergreen Shrub" \
1.317 + ,"Broadleaf Deciduous Temperate Shrub" \
1.318 + ,"Broadleaf Deciduous Boreal Shrub" \
1.319 + ,"C3 Arctic Grass" \
1.320 + ,"C3 Non-Arctic Grass" \
1.321 + ,"C4 Grass" \
1.322 + ,"Corn" \
1.323 +; ,"Wheat" \
1.324 + ,"All Biome" \
1.325 + /)
1.326 + nrow = dimsizes(row_head)
1.327 +
1.328 +; arrays to be passed to table.
1.329 + text4 = new ((/nrow, ncol/),string )
1.330 +
1.331 + do i=0,nrow-2
1.332 + text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)
1.333 + text4(i,1) = sprintf("%.1f",NPP1_g(i))
1.334 + text4(i,2) = sprintf("%.1f",VegC_g(i))
1.335 + text4(i,3) = sprintf("%.1f",LiCwC_g(i))
1.336 + text4(i,4) = sprintf("%.1f",SoilC_g(i))
1.337 + text4(i,5) = sprintf("%.2f",NPP_ratio(i))
1.338 + text4(i,6) = sprintf("%.1f",NEE2_g(i))
1.339 + end do
1.340 +
1.341 +;-------------------------------------------------------
1.342 +; create html table1
1.343 +
1.344 + header_text = "<H1>NEE and Carbon Stocks and Fluxes: Model "+model_name+"</H1>"
1.345 +
1.346 + header = (/"<HTML>" \
1.347 + ,"<HEAD>" \
1.348 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.349 + ,"</HEAD>" \
1.350 + ,header_text \
1.351 + /)
1.352 + footer = "</HTML>"
1.353 +
1.354 + table_header = (/ \
1.355 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.356 + ,"<tr>" \
1.357 + ," <th bgcolor=DDDDDD >Biome Type</th>" \
1.358 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.359 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.360 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.361 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.362 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.363 + ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
1.364 + ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
1.365 + ,"</tr>" \
1.366 + /)
1.367 + table_footer = "</table>"
1.368 + row_header = "<tr>"
1.369 + row_footer = "</tr>"
1.370 +
1.371 + lines = new(50000,string)
1.372 + nline = 0
1.373 +
1.374 + set_line(lines,nline,header)
1.375 + set_line(lines,nline,table_header)
1.376 +
1.377 +;----------------------------
1.378 +;row of table
1.379 +
1.380 + do n = 0,nrow-2
1.381 + set_line(lines,nline,row_header)
1.382 +
1.383 + txt0 = row_head(n)
1.384 + txt1 = text4(n,0)
1.385 + txt2 = text4(n,1)
1.386 + txt3 = text4(n,2)
1.387 + txt4 = text4(n,3)
1.388 + txt5 = text4(n,4)
1.389 + txt6 = text4(n,5)
1.390 + txt7 = text4(n,6)
1.391 +
1.392 + set_line(lines,nline,"<th>"+txt0+"</th>")
1.393 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.394 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.395 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.396 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.397 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.398 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.399 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.400 +
1.401 + set_line(lines,nline,row_footer)
1.402 + end do
1.403 +;----------------------------
1.404 + set_line(lines,nline,table_footer)
1.405 + set_line(lines,nline,footer)
1.406 +
1.407 +; Now write to an HTML file.
1.408 +
1.409 + output_html = "table_carbon_sink1.html"
1.410 +
1.411 + idx = ind(.not.ismissing(lines))
1.412 + if(.not.any(ismissing(idx))) then
1.413 + asciiwrite(output_html,lines(idx))
1.414 + else
1.415 + print ("error?")
1.416 + end if
1.417 +
1.418 + delete (idx)
1.419 +
1.420 + delete (col_head)
1.421 + delete (row_head)
1.422 + delete (text4)
1.423 + delete (table_header)
1.424 +
1.425 +;-----------------------------------------------------------------
1.426 +; html table2 data
1.427 +
1.428 +; column (not including header column)
1.429 +
1.430 + col_head = (/"NPP (PgC/yr)" \
1.431 + ,"VegC (PgC)" \
1.432 + ,"Litter+CWD (PgC)" \
1.433 + ,"SoilC (PgC)" \
1.434 + ,"NEE (PgC/yr)" \
1.435 + ,"NPP timeseries" \
1.436 + ,"NEE timeseries" \
1.437 + ,"Fire timeseries" \
1.438 + /)
1.439 +
1.440 + ncol = dimsizes(col_head)
1.441 +
1.442 +; row (not including header row)
1.443 +
1.444 +; using model biome class:
1.445 + row_head = (/"Not Vegetated" \
1.446 + ,"Needleleaf Evergreen Temperate Tree" \
1.447 + ,"Needleleaf Evergreen Boreal Tree" \
1.448 +; ,"Needleleaf Deciduous Boreal Tree" \
1.449 + ,"Broadleaf Evergreen Tropical Tree" \
1.450 + ,"Broadleaf Evergreen Temperate Tree" \
1.451 + ,"Broadleaf Deciduous Tropical Tree" \
1.452 + ,"Broadleaf Deciduous Temperate Tree" \
1.453 +; ,"Broadleaf Deciduous Boreal Tree" \
1.454 +; ,"Broadleaf Evergreen Shrub" \
1.455 + ,"Broadleaf Deciduous Temperate Shrub" \
1.456 + ,"Broadleaf Deciduous Boreal Shrub" \
1.457 + ,"C3 Arctic Grass" \
1.458 + ,"C3 Non-Arctic Grass" \
1.459 + ,"C4 Grass" \
1.460 + ,"Corn" \
1.461 +; ,"Wheat" \
1.462 + ,"All Biome" \
1.463 + /)
1.464 + nrow = dimsizes(row_head)
1.465 +
1.466 +; arrays to be passed to table.
1.467 + text4 = new ((/nrow, ncol/),string )
1.468 +
1.469 + do i=0,nrow-2
1.470 + text4(i,0) = sprintf("%.1f",NPP1_t(i))
1.471 + text4(i,1) = sprintf("%.1f",VegC_t(i))
1.472 + text4(i,2) = sprintf("%.1f",LiCwC_t(i))
1.473 + text4(i,3) = sprintf("%.1f",SoilC_t(i))
1.474 + text4(i,4) = sprintf("%.1f",NEE2_t(i))
1.475 + text4(i,5) = "<a href=./NPP_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NPP_annual_biome_"+i+".png>annual_plot</a>"
1.476 + text4(i,6) = "<a href=./NEE_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NEE_annual_biome_"+i+".png>annual_plot</a>"
1.477 + text4(i,7) = "<a href=./Fire_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./Fire_annual_biome_"+i+".png>annual_plot</a>"
1.478 + end do
1.479 + text4(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))
1.480 + text4(nrow-1,1) = sprintf("%.1f",sum(VegC_t))
1.481 + text4(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))
1.482 + text4(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))
1.483 + text4(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))
1.484 + text4(nrow-1,5) = "<a href=./NPP_monthly_global.png>monthly_plot</a> <br> <a href=./NPP_annual_global.png>annual_plot</a>"
1.485 + text4(nrow-1,6) = "<a href=./NEE_monthly_global.png>monthly_plot</a> <br> <a href=./NEE_annual_global.png>annual_plot</a>"
1.486 + text4(nrow-1,7) = "<a href=./Fire_monthly_global.png>monthly_plot</a> <br> <a href=./Fire_annual_global.png>annual_plot</a>"
1.487 +
1.488 +;**************************************************
1.489 +; create html table2
1.490 +;**************************************************
1.491 +
1.492 + header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome): Model "+model_name+"</H1>"
1.493 +
1.494 + header = (/"<HTML>" \
1.495 + ,"<HEAD>" \
1.496 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.497 + ,"</HEAD>" \
1.498 + ,header_text \
1.499 + /)
1.500 + footer = "</HTML>"
1.501 +
1.502 + table_header = (/ \
1.503 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.504 + ,"<tr>" \
1.505 + ," <th bgcolor=DDDDDD >Biome Type</th>" \
1.506 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.507 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.508 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.509 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.510 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.511 + ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
1.512 + ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
1.513 + ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
1.514 + ,"</tr>" \
1.515 + /)
1.516 + table_footer = "</table>"
1.517 + row_header = "<tr>"
1.518 + row_footer = "</tr>"
1.519 +
1.520 + lines = new(50000,string)
1.521 + nline = 0
1.522 +
1.523 + set_line(lines,nline,header)
1.524 + set_line(lines,nline,table_header)
1.525 +;-----------------------------------------------
1.526 +;row of table
1.527 +
1.528 + do n = 0,nrow-1
1.529 + set_line(lines,nline,row_header)
1.530 +
1.531 + txt0 = row_head(n)
1.532 + txt1 = text4(n,0)
1.533 + txt2 = text4(n,1)
1.534 + txt3 = text4(n,2)
1.535 + txt4 = text4(n,3)
1.536 + txt5 = text4(n,4)
1.537 + txt6 = text4(n,5)
1.538 + txt7 = text4(n,6)
1.539 + txt8 = text4(n,7)
1.540 +
1.541 + set_line(lines,nline,"<th>"+txt0+"</th>")
1.542 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.543 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.544 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.545 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.546 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.547 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.548 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.549 + set_line(lines,nline,"<th>"+txt8+"</th>")
1.550 +
1.551 + set_line(lines,nline,row_footer)
1.552 + end do
1.553 +;-----------------------------------------------
1.554 + set_line(lines,nline,table_footer)
1.555 + set_line(lines,nline,footer)
1.556 +
1.557 +; Now write to an HTML file.
1.558 +
1.559 + output_html = "table_carbon_sink2.html"
1.560 +
1.561 + idx = ind(.not.ismissing(lines))
1.562 + if(.not.any(ismissing(idx))) then
1.563 + asciiwrite(output_html,lines(idx))
1.564 + else
1.565 + print ("error?")
1.566 + end if
1.567 +
1.568 + delete (idx)
1.569 +
1.570 +;---------------------------------------------------
1.571 +; read data: model, time series
1.572 +
1.573 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.574 + film = model_name2 + "_Fire_C_1979-2004_monthly.nc"
1.575 + fm = addfile (dirm+film,"r")
1.576 +
1.577 + NPP3 = fm->NPP
1.578 + NEE3 = fm->NEE
1.579 + Fire = fm->COL_FIRE_CLOSS
1.580 +
1.581 + delete (fm)
1.582 +
1.583 +; Units for these variables are:
1.584 +
1.585 +;NPP3: g C/m^2/s
1.586 +;NEE3: g C/m^2/s
1.587 +;Fire: g C/m^2/s
1.588 +
1.589 + nsec_per_month = 60*60*24*30
1.590 +
1.591 +; change unit to g C/m^2/month
1.592 +
1.593 + NPP3 = NPP3 * nsec_per_month
1.594 + NEE3 = NEE3 * nsec_per_month
1.595 + Fire = Fire * nsec_per_month
1.596 +
1.597 + data_n = 3
1.598 +
1.599 + dsizes = dimsizes(NPP3)
1.600 + nyear = dsizes(0)
1.601 + nmonth = dsizes(1)
1.602 + ntime = nyear * nmonth
1.603 +
1.604 + year_start = 1979
1.605 + year_end = 2004
1.606 +
1.607 +;*******************************************************************
1.608 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.609 +;********************************************************************
1.610 +
1.611 +; using model biome class
1.612 + nclass = nclass_mod
1.613 +
1.614 + range = fspan(0,nclass,nclass+1)
1.615 +
1.616 +; print (range)
1.617 +; Use this range information to grab all the values in a
1.618 +; particular range, and then take an average.
1.619 +
1.620 + nx = dimsizes(range) - 1
1.621 +
1.622 +;==============================
1.623 +; put data into bins
1.624 +;==============================
1.625 +
1.626 +; using observed biome class
1.627 +; base = ndtooned(classob)
1.628 +; using model biome class
1.629 + base = ndtooned(classmod)
1.630 +
1.631 +; output
1.632 +
1.633 + yvalues = new((/ntime,data_n,nx/),float)
1.634 +
1.635 +; Loop through each range, using base.
1.636 +
1.637 + do i=0,nx-1
1.638 +
1.639 + if (i.ne.(nx-1)) then
1.640 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.641 + else
1.642 + idx = ind(base.ge.range(i))
1.643 + end if
1.644 +
1.645 + do n = 0,data_n-1
1.646 +
1.647 + t = -1
1.648 + do m = 0,nyear-1
1.649 + do k = 0,nmonth-1
1.650 +
1.651 + t = t + 1
1.652 +
1.653 + if (n.eq.0) then
1.654 + data = ndtooned(NPP3(m,k,:,:))
1.655 + end if
1.656 +
1.657 + if (n.eq.1) then
1.658 + data = ndtooned(NEE3(m,k,:,:))
1.659 + end if
1.660 +
1.661 + if (n.eq.2) then
1.662 + data = ndtooned(Fire(m,k,:,:))
1.663 + end if
1.664 +
1.665 +; Calculate average
1.666 +
1.667 + if (.not.any(ismissing(idx))) then
1.668 + yvalues(t,n,i) = sum(data(idx)*area_1d(idx))
1.669 + else
1.670 + yvalues(t,n,i) = yvalues@_FillValue
1.671 + end if
1.672 +
1.673 +;#############################################################
1.674 +; using model biome class:
1.675 +; set the following 4 classes to _FillValue:
1.676 +; (3)Needleleaf Deciduous Boreal Tree,
1.677 +; (8)Broadleaf Deciduous Boreal Tree,
1.678 +; (9)Broadleaf Evergreen Shrub,
1.679 +; (16)Wheat
1.680 +
1.681 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.682 + yvalues(t,n,i) = yvalues@_FillValue
1.683 + end if
1.684 +;#############################################################
1.685 +
1.686 + end do
1.687 + end do
1.688 +
1.689 + delete(data)
1.690 + end do
1.691 +
1.692 + delete(idx)
1.693 + end do
1.694 +
1.695 + delete (base)
1.696 + delete (NPP3)
1.697 + delete (NEE3)
1.698 + delete (Fire)
1.699 +
1.700 +;----------------------------------------------------------------
1.701 +; data for tseries plot
1.702 +
1.703 + yvalues_g = new((/ntime,data_n,n_biome/),float)
1.704 +
1.705 + yvalues_g@units = "TgC/month"
1.706 +
1.707 +; change unit to Tg C/month
1.708 +; change unit from g to Tg (Tera gram)
1.709 + factor_unit = 1.e-12
1.710 +
1.711 + yvalues_g = yvalues(:,:,good) * factor_unit
1.712 +
1.713 +;*******************************************************************
1.714 +; general settings for line plot
1.715 +;*******************************************************************
1.716 +
1.717 +; res
1.718 + res = True
1.719 + res@xyDashPatterns = (/0/) ; make lines solid
1.720 + res@xyLineThicknesses = (/2.0/) ; make lines thicker
1.721 + res@xyLineColors = (/"blue"/) ; line color
1.722 +
1.723 + res@trXMinF = year_start
1.724 + res@trXMaxF = year_end + 1
1.725 +
1.726 + res@vpHeightF = 0.4 ; change aspect ratio of plot
1.727 +; res@vpWidthF = 0.8
1.728 + res@vpWidthF = 0.75
1.729 +
1.730 +; res@gsnMaximize = True
1.731 +
1.732 +;*******************************************************************
1.733 +; (A) 1 component in each biome: monthly
1.734 +;*******************************************************************
1.735 +
1.736 + component = (/"NPP","NEE","Fire"/)
1.737 +
1.738 +; for x-axis in xyplot
1.739 +
1.740 + timeI = new((/ntime/),integer)
1.741 + timeF = new((/ntime/),float)
1.742 + timeI = ispan(1,ntime,1)
1.743 + timeF = year_start + (timeI-1)/12.
1.744 + timeF@long_name = "year"
1.745 +
1.746 + plot_data = new((/ntime/),float)
1.747 + plot_data@long_name = "TgC/month"
1.748 +
1.749 + do n = 0, data_n-1
1.750 + do m = 0, n_biome-1
1.751 +
1.752 + plot_name = component(n)+"_monthly_biome_"+ m
1.753 +
1.754 + wks = gsn_open_wks (plot_type,plot_name)
1.755 +
1.756 + title = component(n)+ ": "+ row_head(m)
1.757 + res@tiMainString = title
1.758 + res@tiMainFontHeightF = 0.025
1.759 +
1.760 + plot_data(:) = yvalues_g(:,n,m)
1.761 +
1.762 + plot=gsn_csm_xy(wks,timeF,plot_data,res)
1.763 +
1.764 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.765 + "rm "+plot_name+"."+plot_type)
1.766 +
1.767 + clear (wks)
1.768 + delete (plot)
1.769 + end do
1.770 + end do
1.771 +
1.772 + do n = 0, data_n-1
1.773 +
1.774 + plot_name = component(n)+"_monthly_global"
1.775 +
1.776 + wks = gsn_open_wks (plot_type,plot_name)
1.777 +
1.778 + title = component(n)+ ": Global"
1.779 + res@tiMainString = title
1.780 + res@tiMainFontHeightF = 0.025
1.781 +
1.782 + do k = 0,ntime-1
1.783 + plot_data(k) = sum(yvalues_g(k,n,:))
1.784 + end do
1.785 +
1.786 + plot=gsn_csm_xy(wks,timeF,plot_data,res)
1.787 +
1.788 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.789 + "rm "+plot_name+"."+plot_type)
1.790 +
1.791 + clear (wks)
1.792 + delete (plot)
1.793 + end do
1.794 +
1.795 + delete (plot_data)
1.796 + delete (timeI)
1.797 + delete (timeF)
1.798 +
1.799 +;*******************************************************************
1.800 +; (B) 1 component in each biome: annually
1.801 +;*******************************************************************
1.802 +
1.803 + yvalues_a = new((/nyear,data_n,n_biome/),float)
1.804 + yvalues_g!0 = "time"
1.805 + yvalues_g!1 = "case"
1.806 + yvalues_g!2 = "record"
1.807 +
1.808 + yvalues_a = month_to_annual(yvalues_g,0)
1.809 +
1.810 + delete (yvalues_g)
1.811 +
1.812 +; for x-axis in xyplot
1.813 +
1.814 + timeI = new((/nyear/),integer)
1.815 + timeF = new((/nyear/),float)
1.816 + timeI = ispan(1,nyear,1)
1.817 + timeF = year_start + (timeI-1)
1.818 + timeF@long_name = "year"
1.819 +
1.820 + plot_data = new((/nyear/),float)
1.821 + plot_data@long_name = "TgC/year"
1.822 +
1.823 + do n = 0, data_n-1
1.824 + do m = 0, n_biome-1
1.825 +
1.826 + plot_name = component(n)+"_annual_biome_"+ m
1.827 +
1.828 + wks = gsn_open_wks (plot_type,plot_name)
1.829 +
1.830 + title = component(n)+ ": "+ row_head(m)
1.831 + res@tiMainString = title
1.832 + res@tiMainFontHeightF = 0.025
1.833 +
1.834 + plot_data(:) = yvalues_a(:,n,m)
1.835 +
1.836 + plot=gsn_csm_xy(wks,timeF,plot_data,res)
1.837 +
1.838 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.839 + "rm "+plot_name+"."+plot_type)
1.840 +
1.841 + clear (wks)
1.842 + delete (plot)
1.843 + end do
1.844 + end do
1.845 +
1.846 + do n = 0, data_n-1
1.847 +
1.848 + plot_name = component(n)+"_annual_global"
1.849 +
1.850 + wks = gsn_open_wks (plot_type,plot_name)
1.851 +
1.852 + title = component(n)+ ": Global"
1.853 + res@tiMainString = title
1.854 + res@tiMainFontHeightF = 0.025
1.855 +
1.856 + do k = 0,nyear-1
1.857 + plot_data(k) = sum(yvalues_a(k,n,:))
1.858 + end do
1.859 +
1.860 + plot=gsn_csm_xy(wks,timeF,plot_data,res)
1.861 +
1.862 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.863 + "rm "+plot_name+"."+plot_type)
1.864 +
1.865 + clear (wks)
1.866 + delete (plot)
1.867 + end do
1.868 +
1.869 +end
1.870 +