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 ;**************************************************************
36 component = (/"Leaf","Wood","Fine_Root","Litter","Coarse_Woody_Debris","Soil"/)
37 n_comp = dimsizes(component)
39 field_pool = (/"LEAFC","WOODC","FROOTC","LITTERC","CWDC","SOILC"/)
40 field_flux = (/"LEAFC_ALLOC","WOODC_ALLOC","FROOTC_ALLOC","LITTERC_LOSS","CWDC_LOSS","SOILC_HR"/)
42 ;************************************************
44 ;************************************************
48 model_name = "i01.06cn"
49 ;model_name = "i01.06casa"
51 ;************************************************
53 ;************************************************
54 ;------------------------------------------------
55 ; read biome data: observed
57 ; biome_name_ob = "MODIS LandCover"
59 ; diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
60 ; filo = "land_class_"+model_grid+".nc"
62 ; fo = addfile(diro+filo,"r")
64 ; classob = tofloat(fo->LAND_CLASS)
68 ; observed data has 20 land-type classes
72 ;---------------------------------------------------
73 ; read biome data: model
75 biome_name_mod = "Model PFT Class"
77 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
78 film = "class_pft_"+model_grid+".nc"
80 fm = addfile(dirm+film,"r")
82 classmod = fm->CLASS_PFT
86 ; model data has 17 land-type classes
90 ;---------------------------------------------------
91 ; read data: model, each component
93 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
94 film = model_name + "_1980-2004_ANN_climo.nc"
96 fm = addfile (dirm+film,"r")
100 pool = fm->$field_pool(k)$
101 flux = fm->$field_flux(k)$
103 ;Units for these variables are:
107 nsec_per_year = 60*60*24*365
109 flux = flux * nsec_per_year
112 ; all the plant pools (leaf, wood, and fine root) and
113 ; coarse woody debris (cwd) and litter pools for
114 ; CASA need to be divided by 1200. The soil flux
115 ; and turnover time are fine and do not need to be adjusted.
117 ; if (k .ne. n_comp-1) then
121 ;*******************************************************************
122 ; Calculate "nice" bins for binning the data in equally spaced ranges
123 ;********************************************************************
125 ; using observed biome class
127 ; using model biome class
131 range = fspan(0,nclassn-1,nclassn)
134 ; Use this range information to grab all the values in a
135 ; particular range, and then take an average.
139 xvalues = new((/2,nx/),float)
140 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
141 dx = xvalues(0,1) - xvalues(0,0) ; range width
142 dx4 = dx/4 ; 1/4 of the range
143 xvalues(1,:) = xvalues(0,:) - dx/5.
145 ;==============================
147 ;==============================
149 ; using observed biome class
150 ; base_1D = ndtooned(classob)
151 ; using model biome class
152 base_1D = ndtooned(classmod)
154 data1_1D = ndtooned(pool)
155 data2_1D = ndtooned(flux)
159 yvalues = new((/2,nx/),float)
160 count = new((/2,nx/),float)
164 ; See if we are doing data1 (nd=0) or data2 (nd=1).
174 ; Loop through each range, using base.
177 if (i.ne.(nr-2)) then
179 ; print("In range ["+range(i)+","+range(i+1)+")")
180 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
183 ; print("In range ["+range(i)+",)")
184 idx = ind(base.ge.range(i))
189 if(.not.any(ismissing(idx))) then
190 yvalues(nd,i) = avg(data(idx))
191 count(nd,i) = dimsizes(idx)
193 yvalues(nd,i) = yvalues@_FillValue
197 ;#############################################################
198 ;using observed biome class:
199 ; set the following 4 classes to _FillValue:
200 ; Water Bodies(0), Urban and Build-Up(13),
201 ; Permenant Snow and Ice(15), Unclassified(17)
203 ; if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
204 ; yvalues(nd,i) = yvalues@_FillValue
207 ;#############################################################
209 ;#############################################################
210 ;using model biome class:
211 ; set the following 4 classes to _FillValue:
212 ; (3)Needleleaf Deciduous Boreal Tree,
213 ; (8)Broadleaf Deciduous Boreal Tree,
214 ; (9)Broadleaf Evergreen Shrub,
217 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
218 yvalues(nd,i) = yvalues@_FillValue
221 ;#############################################################
223 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
225 ; Clean up for next time in loop.
233 ;============================
234 ;compute turnover time
235 ;============================
245 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
252 ;print (dimsizes(uu))
253 ;print (dimsizes(vv))
255 uu_count = u_count(good)
256 vv_count = v_count(good)
258 n_biome = dimsizes(uu)
259 t_biome = new((/n_biome/),float)
263 ;t_biome_avg = avg(t_biome)
264 t_biome_avg = sum(uu*uu_count)/sum(vv*vv_count)
269 ;===========================
270 ; for html table - biome
271 ;===========================
273 output_html = "table_"+component(k)+".html"
275 ; column (not including header column)
277 col_head = (/component(k)+" Flux",component(k)+" Pool",component(k)+" Turnover Time"/)
279 ncol = dimsizes(col_head)
281 ; row (not including header row)
283 ;----------------------------------------------------
284 ; using observed biome class:
285 ; row_head = (/"Evergreen Needleleaf Forests" \
286 ; ,"Evergreen Broadleaf Forests" \
287 ; ,"Deciduous Needleleaf Forest" \
288 ; ,"Deciduous Broadleaf Forests" \
290 ; ,"Closed Bushlands" \
291 ; ,"Open Bushlands" \
292 ; ,"Woody Savannas (S. Hem.)" \
293 ; ,"Savannas (S. Hem.)" \
295 ; ,"Permanent Wetlands" \
297 ; ,"Cropland/Natural Vegetation Mosaic" \
298 ; ,"Barren or Sparsely Vegetated" \
299 ; ,"Woody Savannas (N. Hem.)" \
300 ; ,"Savannas (N. Hem.)" \
304 ;----------------------------------------------------
305 ; using model biome class:
306 row_head = (/"Not Vegetated" \
307 ,"Needleleaf Evergreen Temperate Tree" \
308 ,"Needleleaf Evergreen Boreal Tree" \
309 ; ,"Needleleaf Deciduous Boreal Tree" \
310 ,"Broadleaf Evergreen Tropical Tree" \
311 ,"Broadleaf Evergreen Temperate Tree" \
312 ,"Broadleaf Deciduous Tropical Tree" \
313 ,"Broadleaf Deciduous Temperate Tree" \
314 ; ,"Broadleaf Deciduous Boreal Tree" \
315 ; ,"Broadleaf Evergreen Shrub" \
316 ,"Broadleaf Deciduous Temperate Shrub" \
317 ,"Broadleaf Deciduous Boreal Shrub" \
319 ,"C3 Non-Arctic Grass" \
325 nrow = dimsizes(row_head)
327 ; arrays to be passed to table.
328 text4 = new ((/nrow, ncol/),string )
331 text4(i,0) = sprintf("%.1f",vv(i))
332 text4(i,1) = sprintf("%.1f",uu(i))
333 text4(i,2) = sprintf("%.2f",t_biome(i))
335 text4(nrow-1,0) = "-"
336 text4(nrow-1,1) = "-"
337 text4(nrow-1,2) = sprintf("%.2f",t_biome_avg)
339 ;**************************************************
341 ;**************************************************
343 header_text = "<H1>"+component(k)+" Turnover Time: Model "+model_name+"</H1>"
345 header = (/"<HTML>" \
347 ,"<TITLE>CLAMP metrics</TITLE>" \
354 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
356 ," <th bgcolor=DDDDDD >Biome Class</th>" \
357 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
358 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
359 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
362 table_footer = "</table>"
366 lines = new(50000,string)
369 set_line(lines,nline,header)
370 set_line(lines,nline,table_header)
371 ;-----------------------------------------------
375 set_line(lines,nline,row_header)
382 set_line(lines,nline,"<th>"+txt1+"</th>")
383 set_line(lines,nline,"<th>"+txt2+"</th>")
384 set_line(lines,nline,"<th>"+txt3+"</th>")
385 set_line(lines,nline,"<th>"+txt4+"</th>")
387 set_line(lines,nline,row_footer)
389 ;-----------------------------------------------
390 set_line(lines,nline,table_footer)
391 set_line(lines,nline,footer)
393 ; Now write to an HTML file.
394 idx = ind(.not.ismissing(lines))
395 if(.not.any(ismissing(idx))) then
396 asciiwrite(output_html,lines(idx))