1 ;********************************************************
2 ;using model biome vlass
4 ; required command line input parameters:
5 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
7 ; histogram normalized by rain and compute correleration
8 ;**************************************************************
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
12 ;**************************************************************
13 procedure set_line(lines:string,nline:integer,newlines:string)
15 ; add line to ascci/html file
17 nnewlines = dimsizes(newlines)
18 if(nline+nnewlines-1.ge.dimsizes(lines))
19 print("set_line: bad index, not setting anything.")
22 lines(nline:nline+nnewlines-1) = newlines
23 ; print ("lines = " + lines(nline:nline+nnewlines-1))
24 nline = nline + nnewlines
27 ;**************************************************************
34 ;************************************************
36 ;************************************************
41 model_name1 = "i01.06cn"
42 model_name2 = "i01.10cn"
44 ;---------------------------------------------------
45 ; get biome data: model
47 biome_name_mod = "Model PFT Class"
49 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
50 film = "class_pft_"+model_grid+".nc"
52 fm = addfile(dirm+film,"r")
54 classmod = fm->CLASS_PFT
58 ; model data has 17 land-type classes
62 ;--------------------------------------------------
63 ; get model data: landfrac and area
65 dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/"
67 fm_l = addfile (dirm_l+film_l,"r")
69 landfrac = fm_l->landfrac
72 ; change area from km**2 to m**2
74 ;---------------------------------------------------
75 ; read data: model, group 1
77 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
78 film = model_name1 + "_1980-2004_ANN_climo.nc"
79 fm = addfile (dirm+film,"r")
87 VegC = leafc + woodc + frootc
92 LiCwC = litterc + cwdc
97 ;---------------------------------------------------
98 ; read data: model, group 2
100 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
101 film = model_name2 + "_1990-2004_ANN_climo.nc"
102 fm = addfile (dirm+film,"r")
107 ;---------------------------------------------------
108 ; Units for these variables are:
118 nsec_per_year = 60*60*24*365
120 ; change unit to g C/m^2/year
122 NPP1 = NPP1 * nsec_per_year
123 NPP2 = NPP2 * nsec_per_year
124 NEE2 = NEE2 * nsec_per_year
126 ;---------------------------------------------------
127 ; take into account landfrac
129 area(:,:) = area(:,:) * landfrac(:,:)
130 NPP1(0,:,:) = NPP1(0,:,:) * landfrac(:,:)
131 VegC(0,:,:) = VegC(0,:,:) * landfrac(:,:)
132 LiCwC(0,:,:) = LiCwC(0,:,:) * landfrac(:,:)
133 SoilC(0,:,:) = SoilC(0,:,:) * landfrac(:,:)
134 NPP2(0,:,:) = NPP2(0,:,:) * landfrac(:,:)
135 NEE2(0,:,:) = NEE2(0,:,:) * landfrac(:,:)
139 ;*******************************************************************
140 ; Calculate "nice" bins for binning the data in equally spaced ranges
141 ;********************************************************************
143 ; using model biome class
146 range = fspan(0,nclass,nclass+1)
149 ; Use this range information to grab all the values in a
150 ; particular range, and then take an average.
152 nx = dimsizes(range) - 1
154 ;==============================
156 ;==============================
158 ; using observed biome class
159 ; base = ndtooned(classob)
160 ; using model biome class
161 base = ndtooned(classmod)
165 yvalues = new((/data_n,nx/),float)
166 count = new((/data_n,nx/),float)
171 data = ndtooned(area)
175 data = ndtooned(NPP1)
179 data = ndtooned(VegC)
183 data = ndtooned(LiCwC)
187 data = ndtooned(SoilC)
191 data = ndtooned(NPP2)
195 data = ndtooned(NEE2)
198 ; Loop through each range, using base.
201 if (i.ne.(nx-1)) then
203 ; print("In range ["+range(i)+","+range(i+1)+")")
204 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
207 ; print("In range ["+range(i)+",)")
208 idx = ind(base.ge.range(i))
213 if(.not.any(ismissing(idx))) then
215 yvalues(n,i) = sum(data(idx))
217 yvalues(n,i) = avg(data(idx))
220 count(n,i) = dimsizes(idx)
222 yvalues(n,i) = yvalues@_FillValue
226 ;#############################################################
227 ; using model biome class:
228 ; set the following 4 classes to _FillValue:
229 ; (3)Needleleaf Deciduous Boreal Tree,
230 ; (8)Broadleaf Deciduous Boreal Tree,
231 ; (9)Broadleaf Evergreen Shrub,
234 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
235 yvalues(n,i) = yvalues@_FillValue
238 ;#############################################################
240 ; print(n + ": " + count + " points, avg = " + yvalues(n,i))
257 ;----------------------------------------------------------------
260 good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
284 n_biome = dimsizes(NPP1_g)
286 NPP1_t = new((/n_biome/),float)
287 VegC_t = new((/n_biome/),float)
288 LiCwC_t = new((/n_biome/),float)
289 SoilC_t = new((/n_biome/),float)
290 NEE2_t = new((/n_biome/),float)
291 NPP_ratio = new((/n_biome/),float)
293 NPP_ratio = NPP2_g/NPP1_g
295 ;-----------------------------------------------------------------
298 ; change unit from g to Pg (Peta gram)
301 NPP1_t = NPP1_g * area_g * factor_unit
302 VegC_t = VegC_g * area_g * factor_unit
303 LiCwC_t = LiCwC_g * area_g * factor_unit
304 SoilC_t = SoilC_g * area_g * factor_unit
305 NEE2_t = NEE2_g * area_g * factor_unit
309 ;-------------------------------------------------------------
312 ; column (not including header column)
314 col_head = (/"Area (1.e12m2)" \
317 ,"Litter+CWD (gC/m2)" \
323 ncol = dimsizes(col_head)
325 ; row (not including header row)
327 ; using model biome class:
328 row_head = (/"Not Vegetated" \
329 ,"Needleleaf Evergreen Temperate Tree" \
330 ,"Needleleaf Evergreen Boreal Tree" \
331 ; ,"Needleleaf Deciduous Boreal Tree" \
332 ,"Broadleaf Evergreen Tropical Tree" \
333 ,"Broadleaf Evergreen Temperate Tree" \
334 ,"Broadleaf Deciduous Tropical Tree" \
335 ,"Broadleaf Deciduous Temperate Tree" \
336 ; ,"Broadleaf Deciduous Boreal Tree" \
337 ; ,"Broadleaf Evergreen Shrub" \
338 ,"Broadleaf Deciduous Temperate Shrub" \
339 ,"Broadleaf Deciduous Boreal Shrub" \
341 ,"C3 Non-Arctic Grass" \
347 nrow = dimsizes(row_head)
349 ; arrays to be passed to table.
350 text4 = new ((/nrow, ncol/),string )
353 text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)
354 text4(i,1) = sprintf("%.1f",NPP1_g(i))
355 text4(i,2) = sprintf("%.1f",VegC_g(i))
356 text4(i,3) = sprintf("%.1f",LiCwC_g(i))
357 text4(i,4) = sprintf("%.1f",SoilC_g(i))
358 text4(i,5) = sprintf("%.1f",NPP_ratio(i))
359 text4(i,6) = sprintf("%.1f",NEE2_g(i))
361 text4(nrow-1,0) = "-"
362 text4(nrow-1,1) = "-"
363 text4(nrow-1,2) = "-"
365 ;-------------------------------------------------------
368 header_text = "<H1>NEE and Carbon Stocks and Fluxes: Model "+model_name+"</H1>"
370 header = (/"<HTML>" \
372 ,"<TITLE>CLAMP metrics</TITLE>" \
379 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
381 ," <th bgcolor=DDDDDD >Biome Type</th>" \
382 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
383 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
384 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
385 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
386 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
387 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
388 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
391 table_footer = "</table>"
395 lines = new(50000,string)
398 set_line(lines,nline,header)
399 set_line(lines,nline,table_header)
401 ;----------------------------
405 set_line(lines,nline,row_header)
416 set_line(lines,nline,"<th>"+txt0+"</th>")
417 set_line(lines,nline,"<th>"+txt1+"</th>")
418 set_line(lines,nline,"<th>"+txt2+"</th>")
419 set_line(lines,nline,"<th>"+txt3+"</th>")
420 set_line(lines,nline,"<th>"+txt4+"</th>")
421 set_line(lines,nline,"<th>"+txt5+"</th>")
422 set_line(lines,nline,"<th>"+txt6+"</th>")
423 set_line(lines,nline,"<th>"+txt7+"</th>")
425 set_line(lines,nline,row_footer)
427 ;----------------------------
428 set_line(lines,nline,table_footer)
429 set_line(lines,nline,footer)
431 ; Now write to an HTML file.
433 output_html = "table_carbon_sink1.html"
435 idx = ind(.not.ismissing(lines))
436 if(.not.any(ismissing(idx))) then
437 asciiwrite(output_html,lines(idx))
447 delete (table_header)
449 ;-----------------------------------------------------------------
452 ; column (not including header column)
454 col_head = (/"NPP (PgC/yr)" \
456 ,"Litter+CWD (PgC)" \
464 ncol = dimsizes(col_head)
466 ; row (not including header row)
468 ; using model biome class:
469 row_head = (/"Not Vegetated" \
470 ,"Needleleaf Evergreen Temperate Tree" \
471 ,"Needleleaf Evergreen Boreal Tree" \
472 ; ,"Needleleaf Deciduous Boreal Tree" \
473 ,"Broadleaf Evergreen Tropical Tree" \
474 ,"Broadleaf Evergreen Temperate Tree" \
475 ,"Broadleaf Deciduous Tropical Tree" \
476 ,"Broadleaf Deciduous Temperate Tree" \
477 ; ,"Broadleaf Deciduous Boreal Tree" \
478 ; ,"Broadleaf Evergreen Shrub" \
479 ,"Broadleaf Deciduous Temperate Shrub" \
480 ,"Broadleaf Deciduous Boreal Shrub" \
482 ,"C3 Non-Arctic Grass" \
488 nrow = dimsizes(row_head)
490 ; arrays to be passed to table.
491 text4 = new ((/nrow, ncol/),string )
494 text4(i,0) = sprintf("%.1f",NPP1_t(i))
495 text4(i,1) = sprintf("%.1f",VegC_t(i))
496 text4(i,2) = sprintf("%.1f",LiCwC_t(i))
497 text4(i,3) = sprintf("%.1f",SoilC_t(i))
498 text4(i,4) = sprintf("%.1f",NEE2_t(i))
503 text4(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))
504 text4(nrow-1,1) = sprintf("%.1f",sum(VegC_t))
505 text4(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))
506 text4(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))
507 text4(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))
508 text4(nrow-1,5) = "-"
509 text4(nrow-1,6) = "-"
510 text4(nrow-1,7) = "-"
512 ;**************************************************
514 ;**************************************************
516 header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome): Model "+model_name+"</H1>"
518 header = (/"<HTML>" \
520 ,"<TITLE>CLAMP metrics</TITLE>" \
527 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
529 ," <th bgcolor=DDDDDD >Biome Type</th>" \
530 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
531 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
532 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
533 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
534 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
535 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
536 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
537 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
540 table_footer = "</table>"
544 lines = new(50000,string)
547 set_line(lines,nline,header)
548 set_line(lines,nline,table_header)
549 ;-----------------------------------------------
553 set_line(lines,nline,row_header)
565 set_line(lines,nline,"<th>"+txt0+"</th>")
566 set_line(lines,nline,"<th>"+txt1+"</th>")
567 set_line(lines,nline,"<th>"+txt2+"</th>")
568 set_line(lines,nline,"<th>"+txt3+"</th>")
569 set_line(lines,nline,"<th>"+txt4+"</th>")
570 set_line(lines,nline,"<th>"+txt5+"</th>")
571 set_line(lines,nline,"<th>"+txt6+"</th>")
572 set_line(lines,nline,"<th>"+txt7+"</th>")
573 set_line(lines,nline,"<th>"+txt8+"</th>")
575 set_line(lines,nline,row_footer)
577 ;-----------------------------------------------
578 set_line(lines,nline,table_footer)
579 set_line(lines,nline,footer)
581 ; Now write to an HTML file.
583 output_html = "table_carbon_sink2.html"
585 idx = ind(.not.ismissing(lines))
586 if(.not.any(ismissing(idx))) then
587 asciiwrite(output_html,lines(idx))