Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;********************************************************
2 ; hardwire: flux = flux/1200. (for casa only)
3 ;**************************************************************
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
7 ;**************************************************************
8 procedure set_line(lines:string,nline:integer,newlines:string)
10 ; add line to ascci/html file
12 nnewlines = dimsizes(newlines)
13 if(nline+nnewlines-1.ge.dimsizes(lines))
14 print("set_line: bad index, not setting anything.")
17 lines(nline:nline+nnewlines-1) = newlines
18 ; print ("lines = " + lines(nline:nline+nnewlines-1))
19 nline = nline + nnewlines
22 ;**************************************************************
29 ;------------------------------------------------------
30 ; edit table.html of current model for movel1_vs_model2
32 if (isvar("compare")) then
33 html_name2 = compare+"/table.html"
34 html_new2 = html_name2 +".new"
37 ;------------------------------------------------------
38 ; edit table.html for current model
40 html_name = model_name+"/table.html"
41 html_new = html_name +".new"
43 ;---------------------------------------------------------------
46 component = (/"Leaf","Wood","Fine_Root","Litter","Coarse_Woody_Debris","Soil"/)
47 n_comp = dimsizes(component)
49 field_pool = (/"LEAFC","WOODC","FROOTC","LITTERC","CWDC","SOILC"/)
50 field_flux = (/"LEAFC_ALLOC","WOODC_ALLOC","FROOTC_ALLOC","LITTERC_LOSS","CWDC_LOSS","SOILC_HR"/)
52 ;--------------------------------------------------
55 film_l = "lnd_"+ model_grid +".nc"
56 fm_l = addfile (dirs+film_l,"r")
57 landfrac = fm_l->landfrac
60 ;---------------------------------------------------
61 ; read biome data: model
63 biome_name_mod = "Model PFT Class"
65 film_c = "class_pft_"+ model_grid +".nc"
66 fm_c = addfile (dirs+film_c,"r")
67 classmod = fm_c->CLASS_PFT
71 ; model data has 17 land-type classes
74 ;********************************************************************
75 ; use land-type class to bin the data in equally spaced ranges
76 ;********************************************************************
78 ; using model biome class
81 range = fspan(0,nclass,nclass+1)
83 ; Use this range information to grab all the values in a
84 ; particular range, and then take an average.
86 nx = dimsizes(range) - 1
88 ; for 2 data: pool and flux
91 ; using model biome class
93 base = ndtooned(classmod)
97 yvalues = new((/data_n,nx/),float)
98 count = new((/data_n,nx/),float)
100 ;--------------------------------------------------
101 ; read model data, each component:
103 fm = addfile (dirm+film4,"r")
107 pool = fm->$field_pool(k)$
108 flux = fm->$field_flux(k)$
110 ; Units for these variables are:
114 nsec_per_year = 60*60*24*365
116 flux = flux * nsec_per_year
119 unit_f = "gC/m2/year"
122 ;#############################################################
124 ; all the plant pools (leaf, wood, and fine root) and
125 ; coarse woody debris (cwd) and litter pools for
126 ; CASA need to be divided by 1200. The soil flux
127 ; and turnover time are fine and do not need to be adjusted.
129 if (BGC .eq. "casa") then
130 if (k .ne. n_comp-1) then
134 ;##############################################################
136 ; take into account landfrac
138 pool = pool * conform(pool,landfrac,(/1,2/))
139 flux = flux * conform(flux,landfrac,(/1,2/))
141 ; Loop through each range, using base
145 if (i.ne.(nx-1)) then
146 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
148 idx = ind(base.ge.range(i))
151 ; loop through each dataset
156 data = ndtooned(pool)
160 data = ndtooned(flux)
165 if (.not.any(ismissing(idx))) then
166 yvalues(n,i) = avg(data(idx))
167 count(n,i) = dimsizes(idx)
169 yvalues(n,i) = yvalues@_FillValue
173 ;#############################################################
174 ; using model biome class:
176 ; set the following 4 classes to _FillValue:
177 ; (3)Needleleaf Deciduous Boreal Tree,
178 ; (8)Broadleaf Deciduous Boreal Tree,
179 ; (9)Broadleaf Evergreen Shrub,
182 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
183 yvalues(n,i) = yvalues@_FillValue
186 ;#############################################################
197 ;============================
198 ;compute turnover time
199 ;============================
206 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
210 uu_count = u_count(good)
211 vv_count = v_count(good)
213 n_biome = dimsizes(uu)
214 t_biome = new((/n_biome/),float)
218 t_biome_avg = sum(uu*uu_count)/sum(vv*vv_count)
220 ;===========================
221 ; for html table - biome
222 ;===========================
224 output_html = "table_"+component(k)+".html"
226 ; column (not including header column)
228 col_head = (/component(k)+" Flux",component(k)+" Pool",component(k)+" Turnover Time"/)
230 ncol = dimsizes(col_head)
232 ; row (not including header row)
234 ;----------------------------------------------------
235 ; using model biome class:
236 row_head = (/"Not Vegetated" \
237 ,"Needleleaf Evergreen Temperate Tree" \
238 ,"Needleleaf Evergreen Boreal Tree" \
239 ; ,"Needleleaf Deciduous Boreal Tree" \
240 ,"Broadleaf Evergreen Tropical Tree" \
241 ,"Broadleaf Evergreen Temperate Tree" \
242 ,"Broadleaf Deciduous Tropical Tree" \
243 ,"Broadleaf Deciduous Temperate Tree" \
244 ; ,"Broadleaf Deciduous Boreal Tree" \
245 ; ,"Broadleaf Evergreen Shrub" \
246 ,"Broadleaf Deciduous Temperate Shrub" \
247 ,"Broadleaf Deciduous Boreal Shrub" \
249 ,"C3 Non-Arctic Grass" \
255 nrow = dimsizes(row_head)
257 ; arrays to be passed to table.
258 text = new ((/nrow, ncol/),string )
261 text(i,0) = sprintf("%.1f",vv(i))
262 text(i,1) = sprintf("%.1f",uu(i))
263 text(i,2) = sprintf("%.2f",t_biome(i))
267 text(nrow-1,2) = sprintf("%.2f",t_biome_avg)
269 ;**************************************************
271 ;**************************************************
273 header_text = "<H1>"+component(k)+" Turnover Time: Model "+model_name+"</H1>"
275 header = (/"<HTML>" \
277 ,"<TITLE>CLAMP metrics</TITLE>" \
284 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
286 ," <th bgcolor=DDDDDD >Biome Class</th>" \
287 ," <th bgcolor=DDDDDD >"+col_head(0)+"<br>("+unit_f+")</th>" \
288 ," <th bgcolor=DDDDDD >"+col_head(1)+"<br>("+unit_p+")</th>" \
289 ," <th bgcolor=DDDDDD >"+col_head(2)+"<br>("+unit_t+")</th>" \
292 table_footer = "</table>"
296 lines = new(50000,string)
299 set_line(lines,nline,header)
300 set_line(lines,nline,table_header)
301 ;-----------------------------------------------
305 set_line(lines,nline,row_header)
312 set_line(lines,nline,"<th>"+txt1+"</th>")
313 set_line(lines,nline,"<th>"+txt2+"</th>")
314 set_line(lines,nline,"<th>"+txt3+"</th>")
315 set_line(lines,nline,"<th>"+txt4+"</th>")
317 set_line(lines,nline,row_footer)
319 ;-----------------------------------------------
320 set_line(lines,nline,table_footer)
321 set_line(lines,nline,footer)
323 ; Now write to an HTML file
325 idx = ind(.not.ismissing(lines))
326 if(.not.any(ismissing(idx))) then
327 asciiwrite(output_html,lines(idx))
342 ;***************************************************************************
343 ; output plot and html
344 ;***************************************************************************
345 output_dir = model_name+"/turnover"
347 system("mv *.html " + output_dir)
348 ;******************************