1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/carbon_sink/11.table1.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,447 @@
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 +; read data: model, group 1
1.79 +
1.80 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.81 + film = model_name1 + "_1980-2004_ANN_climo.nc"
1.82 + fm = addfile (dirm+film,"r")
1.83 +
1.84 + NPP1 = fm->NPP
1.85 +
1.86 + leafc = fm->LEAFC
1.87 + woodc = fm->WOODC
1.88 + frootc = fm->FROOTC
1.89 + VegC = leafc
1.90 + VegC = leafc + woodc + frootc
1.91 +
1.92 + litterc = fm->LITTERC
1.93 + cwdc = fm->CWDC
1.94 + LiCwC = litterc
1.95 + LiCwC = litterc + cwdc
1.96 +
1.97 + SoilC = fm->SOILC
1.98 +
1.99 + delete (fm)
1.100 +;---------------------------------------------------
1.101 +; read data: model, group 2
1.102 +
1.103 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.104 + film = model_name2 + "_1990-2004_ANN_climo.nc"
1.105 + fm = addfile (dirm+film,"r")
1.106 +
1.107 + NPP2 = fm->NPP
1.108 + NEE2 = fm->NEE
1.109 +
1.110 +;---------------------------------------------------
1.111 +; Units for these variables are:
1.112 +
1.113 +;NPP1: g C/m^2/s
1.114 +;NPP2: g C/m^2/s
1.115 +;NEE2: g C/m^2/s
1.116 +
1.117 +;VegC: g C/m^2
1.118 +;LiCwC: g C/m^2
1.119 +;SoilC: g C/m^2
1.120 +
1.121 + nsec_per_year = 60*60*24*365
1.122 +
1.123 +; change unit to g C/m^2/year
1.124 +
1.125 + NPP1 = NPP1 * nsec_per_year
1.126 + NPP2 = NPP2 * nsec_per_year
1.127 + NEE2 = NEE2 * nsec_per_year
1.128 +
1.129 +;---------------------------------------------------
1.130 +; take into account landfrac
1.131 +
1.132 + area(:,:) = area(:,:) * landfrac(:,:)
1.133 + NPP1(0,:,:) = NPP1(0,:,:) * landfrac(:,:)
1.134 + VegC(0,:,:) = VegC(0,:,:) * landfrac(:,:)
1.135 + LiCwC(0,:,:) = LiCwC(0,:,:) * landfrac(:,:)
1.136 + SoilC(0,:,:) = SoilC(0,:,:) * landfrac(:,:)
1.137 + NPP2(0,:,:) = NPP2(0,:,:) * landfrac(:,:)
1.138 + NEE2(0,:,:) = NEE2(0,:,:) * landfrac(:,:)
1.139 +
1.140 + data_n = 7
1.141 +
1.142 +;*******************************************************************
1.143 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.144 +;********************************************************************
1.145 +
1.146 +; using model biome class
1.147 + nclass = nclass_mod
1.148 +
1.149 + range = fspan(0,nclass,nclass+1)
1.150 +
1.151 +; print (range)
1.152 +; Use this range information to grab all the values in a
1.153 +; particular range, and then take an average.
1.154 +
1.155 + nx = dimsizes(range) - 1
1.156 +
1.157 +;==============================
1.158 +; put data into bins
1.159 +;==============================
1.160 +
1.161 +; using observed biome class
1.162 +; base = ndtooned(classob)
1.163 +; using model biome class
1.164 + base = ndtooned(classmod)
1.165 +
1.166 +; output
1.167 +
1.168 + yvalues = new((/data_n,nx/),float)
1.169 + count = new((/data_n,nx/),float)
1.170 +
1.171 + do n = 0,data_n-1
1.172 +
1.173 + if(n.eq.0) then
1.174 + data = ndtooned(area)
1.175 + end if
1.176 +
1.177 + if(n.eq.1) then
1.178 + data = ndtooned(NPP1)
1.179 + end if
1.180 +
1.181 + if(n.eq.2) then
1.182 + data = ndtooned(VegC)
1.183 + end if
1.184 +
1.185 + if(n.eq.3) then
1.186 + data = ndtooned(LiCwC)
1.187 + end if
1.188 +
1.189 + if(n.eq.4) then
1.190 + data = ndtooned(SoilC)
1.191 + end if
1.192 +
1.193 + if(n.eq.5) then
1.194 + data = ndtooned(NPP2)
1.195 + end if
1.196 +
1.197 + if(n.eq.6) then
1.198 + data = ndtooned(NEE2)
1.199 + end if
1.200 +
1.201 +; Loop through each range, using base.
1.202 +
1.203 + do i=0,nx-1
1.204 + if (i.ne.(nx-1)) then
1.205 +; print("")
1.206 +; print("In range ["+range(i)+","+range(i+1)+")")
1.207 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.208 + else
1.209 +; print("")
1.210 +; print("In range ["+range(i)+",)")
1.211 + idx = ind(base.ge.range(i))
1.212 + end if
1.213 +
1.214 +; Calculate average
1.215 +
1.216 + if(.not.any(ismissing(idx))) then
1.217 + if (n.eq.0) then
1.218 + yvalues(n,i) = sum(data(idx))
1.219 + else
1.220 + yvalues(n,i) = avg(data(idx))
1.221 + end if
1.222 +
1.223 + count(n,i) = dimsizes(idx)
1.224 + else
1.225 + yvalues(n,i) = yvalues@_FillValue
1.226 + count(n,i) = 0
1.227 + end if
1.228 +
1.229 +;#############################################################
1.230 +; using model biome class:
1.231 +; set the following 4 classes to _FillValue:
1.232 +; (3)Needleleaf Deciduous Boreal Tree,
1.233 +; (8)Broadleaf Deciduous Boreal Tree,
1.234 +; (9)Broadleaf Evergreen Shrub,
1.235 +; (16)Wheat
1.236 +
1.237 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.238 + yvalues(n,i) = yvalues@_FillValue
1.239 + count(n,i) = 0
1.240 + end if
1.241 +;#############################################################
1.242 +
1.243 +; print(n + ": " + count + " points, avg = " + yvalues(n,i))
1.244 +
1.245 + delete(idx)
1.246 + end do
1.247 +
1.248 + delete(data)
1.249 + end do
1.250 +
1.251 + delete (base)
1.252 + delete (area)
1.253 + delete (NPP1)
1.254 + delete (VegC)
1.255 + delete (LiCwC)
1.256 + delete (SoilC)
1.257 + delete (NPP2)
1.258 + delete (NEE2)
1.259 +
1.260 +;----------------------------------------------------------------
1.261 +; data for table1
1.262 +
1.263 + good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
1.264 +;print (good)
1.265 +
1.266 + w = yvalues(0,:)
1.267 + area_g = w(good)
1.268 +
1.269 + w = yvalues(1,:)
1.270 + NPP1_g = w(good)
1.271 +
1.272 + w = yvalues(2,:)
1.273 + VegC_g = w(good)
1.274 +
1.275 + w = yvalues(3,:)
1.276 + LiCwC_g = w(good)
1.277 +
1.278 + w = yvalues(4,:)
1.279 + SoilC_g = w(good)
1.280 +
1.281 + w = yvalues(5,:)
1.282 + NPP2_g = w(good)
1.283 +
1.284 + w = yvalues(6,:)
1.285 + NEE2_g = w(good)
1.286 +
1.287 + n_biome = dimsizes(NPP1_g)
1.288 +
1.289 + NPP1_t = new((/n_biome/),float)
1.290 + VegC_t = new((/n_biome/),float)
1.291 + LiCwC_t = new((/n_biome/),float)
1.292 + SoilC_t = new((/n_biome/),float)
1.293 + NEE2_t = new((/n_biome/),float)
1.294 + NPP_ratio = new((/n_biome/),float)
1.295 +
1.296 + NPP_ratio = NPP2_g/NPP1_g
1.297 +
1.298 +;-----------------------------------------------------------------
1.299 +; data for table2
1.300 +
1.301 +; change unit from g to Pg (Peta gram)
1.302 + factor_unit = 1.e-15
1.303 +
1.304 + NPP1_t = NPP1_g * area_g * factor_unit
1.305 + VegC_t = VegC_g * area_g * factor_unit
1.306 + LiCwC_t = LiCwC_g * area_g * factor_unit
1.307 + SoilC_t = SoilC_g * area_g * factor_unit
1.308 + NEE2_t = NEE2_g * area_g * factor_unit
1.309 +
1.310 + print (NPP1_t)
1.311 +
1.312 +;-------------------------------------------------------------
1.313 +; html table1 data
1.314 +
1.315 +; column (not including header column)
1.316 +
1.317 + col_head = (/"Area (1.e12m2)" \
1.318 + ,"NPP (gC/m2/yr)" \
1.319 + ,"VegC (gC/m2)" \
1.320 + ,"Litter+CWD (gC/m2)" \
1.321 + ,"SoilC (gC/m2)" \
1.322 + ,"NPP ratio" \
1.323 + ,"NEE (gC/m2/yr)" \
1.324 + /)
1.325 +
1.326 + ncol = dimsizes(col_head)
1.327 +
1.328 +; row (not including header row)
1.329 +
1.330 +; using model biome class:
1.331 + row_head = (/"Not Vegetated" \
1.332 + ,"Needleleaf Evergreen Temperate Tree" \
1.333 + ,"Needleleaf Evergreen Boreal Tree" \
1.334 +; ,"Needleleaf Deciduous Boreal Tree" \
1.335 + ,"Broadleaf Evergreen Tropical Tree" \
1.336 + ,"Broadleaf Evergreen Temperate Tree" \
1.337 + ,"Broadleaf Deciduous Tropical Tree" \
1.338 + ,"Broadleaf Deciduous Temperate Tree" \
1.339 +; ,"Broadleaf Deciduous Boreal Tree" \
1.340 +; ,"Broadleaf Evergreen Shrub" \
1.341 + ,"Broadleaf Deciduous Temperate Shrub" \
1.342 + ,"Broadleaf Deciduous Boreal Shrub" \
1.343 + ,"C3 Arctic Grass" \
1.344 + ,"C3 Non-Arctic Grass" \
1.345 + ,"C4 Grass" \
1.346 + ,"Corn" \
1.347 +; ,"Wheat" \
1.348 + ,"All Biome" \
1.349 + /)
1.350 + nrow = dimsizes(row_head)
1.351 +
1.352 +; arrays to be passed to table.
1.353 + text4 = new ((/nrow, ncol/),string )
1.354 +
1.355 + do i=0,nrow-2
1.356 + text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)
1.357 + text4(i,1) = sprintf("%.1f",NPP1_g(i))
1.358 + text4(i,2) = sprintf("%.1f",VegC_g(i))
1.359 + text4(i,3) = sprintf("%.1f",LiCwC_g(i))
1.360 + text4(i,4) = sprintf("%.1f",SoilC_g(i))
1.361 + text4(i,5) = sprintf("%.1f",NPP_ratio(i))
1.362 + text4(i,6) = sprintf("%.1f",NEE2_g(i))
1.363 + end do
1.364 + text4(nrow-1,0) = "-"
1.365 + text4(nrow-1,1) = "-"
1.366 + text4(nrow-1,2) = "-"
1.367 +
1.368 +;-------------------------------------------------------
1.369 +; create html table1
1.370 +
1.371 + header_text = "<H1>NEE and Carbon Stocks and Fluxes: Model "+model_name+"</H1>"
1.372 +
1.373 + header = (/"<HTML>" \
1.374 + ,"<HEAD>" \
1.375 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.376 + ,"</HEAD>" \
1.377 + ,header_text \
1.378 + /)
1.379 + footer = "</HTML>"
1.380 +
1.381 + table_header = (/ \
1.382 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.383 + ,"<tr>" \
1.384 + ," <th bgcolor=DDDDDD >Biome Type</th>" \
1.385 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.386 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.387 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.388 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.389 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.390 + ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
1.391 + ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
1.392 + ,"</tr>" \
1.393 + /)
1.394 + table_footer = "</table>"
1.395 + row_header = "<tr>"
1.396 + row_footer = "</tr>"
1.397 +
1.398 + lines = new(50000,string)
1.399 + nline = 0
1.400 +
1.401 + set_line(lines,nline,header)
1.402 + set_line(lines,nline,table_header)
1.403 +
1.404 +;----------------------------
1.405 +;row of table
1.406 +
1.407 + do n = 0,nrow-2
1.408 + set_line(lines,nline,row_header)
1.409 +
1.410 + txt0 = row_head(n)
1.411 + txt1 = text4(n,0)
1.412 + txt2 = text4(n,1)
1.413 + txt3 = text4(n,2)
1.414 + txt4 = text4(n,3)
1.415 + txt5 = text4(n,4)
1.416 + txt6 = text4(n,5)
1.417 + txt7 = text4(n,6)
1.418 +
1.419 + set_line(lines,nline,"<th>"+txt0+"</th>")
1.420 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.421 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.422 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.423 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.424 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.425 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.426 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.427 +
1.428 + set_line(lines,nline,row_footer)
1.429 + end do
1.430 +;----------------------------
1.431 + set_line(lines,nline,table_footer)
1.432 + set_line(lines,nline,footer)
1.433 +
1.434 +; Now write to an HTML file.
1.435 +
1.436 + output_html = "table_carbon_sink1.html"
1.437 +
1.438 + idx = ind(.not.ismissing(lines))
1.439 + if(.not.any(ismissing(idx))) then
1.440 + asciiwrite(output_html,lines(idx))
1.441 + else
1.442 + print ("error?")
1.443 + end if
1.444 +
1.445 + delete (col_head)
1.446 + delete (row_head)
1.447 + delete (text4)
1.448 +
1.449 +end
1.450 +