1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/turnover/03.biome_model.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,404 @@
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 +;components
1.38 +
1.39 + component = (/"Leaf","Wood","Fine_Root","Litter","Coarse_Woody_Debris","Soil"/)
1.40 + n_comp = dimsizes(component)
1.41 +
1.42 + field_pool = (/"LEAFC","WOODC","FROOTC","LITTERC","CWDC","SOILC"/)
1.43 + field_flux = (/"LEAFC_ALLOC","WOODC_ALLOC","FROOTC_ALLOC","LITTERC_LOSS","CWDC_LOSS","SOILC_HR"/)
1.44 +
1.45 +;************************************************
1.46 +; model name and grid
1.47 +;************************************************
1.48 +
1.49 + model_grid = "T42"
1.50 +
1.51 + model_name = "i01.06cn"
1.52 +;model_name = "i01.06casa"
1.53 +
1.54 +;************************************************
1.55 +; read data: ob
1.56 +;************************************************
1.57 +;------------------------------------------------
1.58 +; read biome data: observed
1.59 +
1.60 +; biome_name_ob = "MODIS LandCover"
1.61 +
1.62 +; diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
1.63 +; filo = "land_class_"+model_grid+".nc"
1.64 +
1.65 +; fo = addfile(diro+filo,"r")
1.66 +
1.67 +; classob = tofloat(fo->LAND_CLASS)
1.68 +
1.69 +; delete (fo)
1.70 +
1.71 +; observed data has 20 land-type classes
1.72 +
1.73 +; nclass_ob = 20
1.74 +
1.75 +;---------------------------------------------------
1.76 +; read biome data: model
1.77 +
1.78 + biome_name_mod = "Model PFT Class"
1.79 +
1.80 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.81 + film = "class_pft_"+model_grid+".nc"
1.82 +
1.83 + fm = addfile(dirm+film,"r")
1.84 +
1.85 + classmod = fm->CLASS_PFT
1.86 +
1.87 + delete (fm)
1.88 +
1.89 +; model data has 17 land-type classes
1.90 +
1.91 + nclass_mod = 17
1.92 +
1.93 +;---------------------------------------------------
1.94 +; read data: model, each component
1.95 +
1.96 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.97 + film = model_name + "_1980-2004_ANN_climo.nc"
1.98 +
1.99 + fm = addfile (dirm+film,"r")
1.100 +
1.101 + do k = 0,n_comp-1
1.102 +
1.103 + pool = fm->$field_pool(k)$
1.104 + flux = fm->$field_flux(k)$
1.105 +
1.106 +;Units for these variables are:
1.107 +;pool: g C/m^2
1.108 +;flux: g C/m^2/s
1.109 +
1.110 + nsec_per_year = 60*60*24*365
1.111 +
1.112 + flux = flux * nsec_per_year
1.113 +
1.114 +; casa only
1.115 +; all the plant pools (leaf, wood, and fine root) and
1.116 +; coarse woody debris (cwd) and litter pools for
1.117 +; CASA need to be divided by 1200. The soil flux
1.118 +; and turnover time are fine and do not need to be adjusted.
1.119 +
1.120 +; if (k .ne. n_comp-1) then
1.121 +; flux = flux/1200.
1.122 +; end if
1.123 +
1.124 +;*******************************************************************
1.125 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.126 +;********************************************************************
1.127 +
1.128 +; using observed biome class
1.129 +; nclass = nclass_ob
1.130 +; using model biome class
1.131 + nclass = nclass_mod
1.132 +
1.133 + nclassn = nclass + 1
1.134 + range = fspan(0,nclassn-1,nclassn)
1.135 +; print (range)
1.136 +
1.137 +; Use this range information to grab all the values in a
1.138 +; particular range, and then take an average.
1.139 +
1.140 + nr = dimsizes(range)
1.141 + nx = nr-1
1.142 + xvalues = new((/2,nx/),float)
1.143 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.144 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.145 + dx4 = dx/4 ; 1/4 of the range
1.146 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.147 +
1.148 +;==============================
1.149 +; put data into bins
1.150 +;==============================
1.151 +
1.152 +; using observed biome class
1.153 +; base_1D = ndtooned(classob)
1.154 +; using model biome class
1.155 + base_1D = ndtooned(classmod)
1.156 +
1.157 + data1_1D = ndtooned(pool)
1.158 + data2_1D = ndtooned(flux)
1.159 +
1.160 +; output
1.161 +
1.162 + yvalues = new((/2,nx/),float)
1.163 + count = new((/2,nx/),float)
1.164 +
1.165 + do nd=0,1
1.166 +
1.167 +; See if we are doing data1 (nd=0) or data2 (nd=1).
1.168 +
1.169 + base = base_1D
1.170 +
1.171 + if(nd.eq.0) then
1.172 + data = data1_1D
1.173 + else
1.174 + data = data2_1D
1.175 + end if
1.176 +
1.177 +; Loop through each range, using base.
1.178 +
1.179 + do i=0,nr-2
1.180 + if (i.ne.(nr-2)) then
1.181 +; print("")
1.182 +; print("In range ["+range(i)+","+range(i+1)+")")
1.183 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.184 + else
1.185 +; print("")
1.186 +; print("In range ["+range(i)+",)")
1.187 + idx = ind(base.ge.range(i))
1.188 + end if
1.189 +
1.190 +; Calculate average
1.191 +
1.192 + if(.not.any(ismissing(idx))) then
1.193 + yvalues(nd,i) = avg(data(idx))
1.194 + count(nd,i) = dimsizes(idx)
1.195 + else
1.196 + yvalues(nd,i) = yvalues@_FillValue
1.197 + count(nd,i) = 0
1.198 + end if
1.199 +
1.200 +;#############################################################
1.201 +;using observed biome class:
1.202 +; set the following 4 classes to _FillValue:
1.203 +; Water Bodies(0), Urban and Build-Up(13),
1.204 +; Permenant Snow and Ice(15), Unclassified(17)
1.205 +
1.206 +; if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
1.207 +; yvalues(nd,i) = yvalues@_FillValue
1.208 +; count(nd,i) = 0
1.209 +; end if
1.210 +;#############################################################
1.211 +
1.212 +;#############################################################
1.213 +;using model biome class:
1.214 +; set the following 4 classes to _FillValue:
1.215 +; (3)Needleleaf Deciduous Boreal Tree,
1.216 +; (8)Broadleaf Deciduous Boreal Tree,
1.217 +; (9)Broadleaf Evergreen Shrub,
1.218 +; (16)Wheat
1.219 +
1.220 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.221 + yvalues(nd,i) = yvalues@_FillValue
1.222 + count(nd,i) = 0
1.223 + end if
1.224 +;#############################################################
1.225 +
1.226 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.227 +
1.228 +; Clean up for next time in loop.
1.229 +
1.230 + delete(idx)
1.231 + end do
1.232 +
1.233 + delete(data)
1.234 + end do
1.235 +
1.236 +;============================
1.237 +;compute turnover time
1.238 +;============================
1.239 +
1.240 + u = yvalues(0,:)
1.241 + v = yvalues(1,:)
1.242 + u_count = count(0,:)
1.243 + v_count = count(1,:)
1.244 +
1.245 +;print (dimsizes(u))
1.246 +;print (dimsizes(v))
1.247 +
1.248 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.249 +
1.250 +;print (good)
1.251 +
1.252 + uu = u(good)
1.253 + vv = v(good)
1.254 +
1.255 +;print (dimsizes(uu))
1.256 +;print (dimsizes(vv))
1.257 +
1.258 + uu_count = u_count(good)
1.259 + vv_count = v_count(good)
1.260 +
1.261 + n_biome = dimsizes(uu)
1.262 + t_biome = new((/n_biome/),float)
1.263 +
1.264 + t_biome = uu/vv
1.265 +
1.266 +;t_biome_avg = avg(t_biome)
1.267 + t_biome_avg = sum(uu*uu_count)/sum(vv*vv_count)
1.268 +
1.269 +;print (t_biome)
1.270 +;print (t_biome_avg)
1.271 +
1.272 +;===========================
1.273 +; for html table - biome
1.274 +;===========================
1.275 +
1.276 + output_html = "table_"+component(k)+".html"
1.277 +
1.278 +; column (not including header column)
1.279 +
1.280 + col_head = (/component(k)+" Flux",component(k)+" Pool",component(k)+" Turnover Time"/)
1.281 +
1.282 + ncol = dimsizes(col_head)
1.283 +
1.284 +; row (not including header row)
1.285 +
1.286 +;----------------------------------------------------
1.287 +; using observed biome class:
1.288 +; row_head = (/"Evergreen Needleleaf Forests" \
1.289 +; ,"Evergreen Broadleaf Forests" \
1.290 +; ,"Deciduous Needleleaf Forest" \
1.291 +; ,"Deciduous Broadleaf Forests" \
1.292 +; ,"Mixed Forests" \
1.293 +; ,"Closed Bushlands" \
1.294 +; ,"Open Bushlands" \
1.295 +; ,"Woody Savannas (S. Hem.)" \
1.296 +; ,"Savannas (S. Hem.)" \
1.297 +; ,"Grasslands" \
1.298 +; ,"Permanent Wetlands" \
1.299 +; ,"Croplands" \
1.300 +; ,"Cropland/Natural Vegetation Mosaic" \
1.301 +; ,"Barren or Sparsely Vegetated" \
1.302 +; ,"Woody Savannas (N. Hem.)" \
1.303 +; ,"Savannas (N. Hem.)" \
1.304 +; ,"All Biome" \
1.305 +; /)
1.306 +
1.307 +;----------------------------------------------------
1.308 +; using model biome class:
1.309 + row_head = (/"Not Vegetated" \
1.310 + ,"Needleleaf Evergreen Temperate Tree" \
1.311 + ,"Needleleaf Evergreen Boreal Tree" \
1.312 +; ,"Needleleaf Deciduous Boreal Tree" \
1.313 + ,"Broadleaf Evergreen Tropical Tree" \
1.314 + ,"Broadleaf Evergreen Temperate Tree" \
1.315 + ,"Broadleaf Deciduous Tropical Tree" \
1.316 + ,"Broadleaf Deciduous Temperate Tree" \
1.317 +; ,"Broadleaf Deciduous Boreal Tree" \
1.318 +; ,"Broadleaf Evergreen Shrub" \
1.319 + ,"Broadleaf Deciduous Temperate Shrub" \
1.320 + ,"Broadleaf Deciduous Boreal Shrub" \
1.321 + ,"C3 Arctic Grass" \
1.322 + ,"C3 Non-Arctic Grass" \
1.323 + ,"C4 Grass" \
1.324 + ,"Corn" \
1.325 +; ,"Wheat" \
1.326 + ,"All Biome" \
1.327 + /)
1.328 + nrow = dimsizes(row_head)
1.329 +
1.330 +; arrays to be passed to table.
1.331 + text4 = new ((/nrow, ncol/),string )
1.332 +
1.333 + do i=0,nrow-2
1.334 + text4(i,0) = sprintf("%.1f",vv(i))
1.335 + text4(i,1) = sprintf("%.1f",uu(i))
1.336 + text4(i,2) = sprintf("%.2f",t_biome(i))
1.337 + end do
1.338 + text4(nrow-1,0) = "-"
1.339 + text4(nrow-1,1) = "-"
1.340 + text4(nrow-1,2) = sprintf("%.2f",t_biome_avg)
1.341 +
1.342 +;**************************************************
1.343 +; html table
1.344 +;**************************************************
1.345 +
1.346 + header_text = "<H1>"+component(k)+" Turnover Time: Model "+model_name+"</H1>"
1.347 +
1.348 + header = (/"<HTML>" \
1.349 + ,"<HEAD>" \
1.350 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.351 + ,"</HEAD>" \
1.352 + ,header_text \
1.353 + /)
1.354 + footer = "</HTML>"
1.355 +
1.356 + table_header = (/ \
1.357 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.358 + ,"<tr>" \
1.359 + ," <th bgcolor=DDDDDD >Biome Class</th>" \
1.360 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.361 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.362 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.363 + ,"</tr>" \
1.364 + /)
1.365 + table_footer = "</table>"
1.366 + row_header = "<tr>"
1.367 + row_footer = "</tr>"
1.368 +
1.369 + lines = new(50000,string)
1.370 + nline = 0
1.371 +
1.372 + set_line(lines,nline,header)
1.373 + set_line(lines,nline,table_header)
1.374 +;-----------------------------------------------
1.375 +;row of table
1.376 +
1.377 + do n = 0,nrow-1
1.378 + set_line(lines,nline,row_header)
1.379 +
1.380 + txt1 = row_head(n)
1.381 + txt2 = text4(n,0)
1.382 + txt3 = text4(n,1)
1.383 + txt4 = text4(n,2)
1.384 +
1.385 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.386 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.387 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.388 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.389 +
1.390 + set_line(lines,nline,row_footer)
1.391 + end do
1.392 +;-----------------------------------------------
1.393 + set_line(lines,nline,table_footer)
1.394 + set_line(lines,nline,footer)
1.395 +
1.396 +; Now write to an HTML file.
1.397 + idx = ind(.not.ismissing(lines))
1.398 + if(.not.any(ismissing(idx))) then
1.399 + asciiwrite(output_html,lines(idx))
1.400 + else
1.401 + print ("error?")
1.402 + end if
1.403 +
1.404 + delete (idx)
1.405 +end do
1.406 +end
1.407 +