1 ;********************************************************
4 ; required command line input parameters:
5 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
7 ;**************************************************************
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
11 ;**************************************************************
12 procedure set_line(lines:string,nline:integer,newlines:string)
14 ; add line to ascci/html file
16 nnewlines = dimsizes(newlines)
17 if(nline+nnewlines-1.ge.dimsizes(lines))
18 print("set_line: bad index, not setting anything.")
21 lines(nline:nline+nnewlines-1) = newlines
22 ; print ("lines = " + lines(nline:nline+nnewlines-1))
23 nline = nline + nnewlines
26 ;**************************************************************
33 ;************************************************
35 ;************************************************
41 ;model_name_i = "i01.07cn"
42 ;model_name_f = "i01.10cn"
44 model_name_i = "i01.07casa"
45 model_name_f = "i01.10casa"
47 model_name = model_name_f
49 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
50 film_i = model_name_i + "_1990-2004_ANN_climo.nc"
51 film_f = model_name_f + "_1990-2004_ANN_climo.nc"
53 fm_i = addfile (dirm+film_i,"r")
54 fm_f = addfile (dirm+film_f,"r")
65 ;Units for these variables are:
68 nsec_per_year = 60*60*24*365
70 npp_i = npp_i * nsec_per_year
71 npp_f = npp_f * nsec_per_year
73 ;===================================================
74 ; read data: observed at stations
75 ;===================================================
77 station = (/"DukeFACE" \
83 lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
84 lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
85 lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
88 n_sta = dimsizes(station)
89 beta_4_ob = new((/n_sta/),float)
92 ;===================================================
93 ; get model data at station
94 ;===================================================
96 npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
98 npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
103 ;============================
105 ;============================
107 beta_4 = new((/n_sta/),float)
109 beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
111 beta_4_avg = avg(beta_4)
116 ;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.
118 bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
119 M_beta = (1. - (bias/n_sta))*3.
123 ;=========================
124 ; for html table - station
125 ;=========================
127 output_html = "table_station.html"
129 ; column (not including header column)
131 col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
133 ncol = dimsizes(col_head)
135 ; row (not including header row)
136 row_head = (/"DukeFACE" \
142 nrow = dimsizes(row_head)
144 ; arrays to be passed to table.
145 text4 = new ((/nrow, ncol/),string )
148 text4(i,0) = sprintf("%.1f",lat_ob(i))
149 text4(i,1) = sprintf("%.1f",lon_ob(i))
150 text4(i,2) = sprintf("%.1f",co2_i)
151 text4(i,3) = sprintf("%.1f",co2_f)
152 text4(i,4) = sprintf("%.1f",npp_i_4(0,i))
153 text4(i,5) = sprintf("%.1f",npp_f_4(0,i))
154 text4(i,6) = sprintf("%.2f",beta_4(i))
157 text4(nrow-1,0) = "-"
158 text4(nrow-1,1) = "-"
159 text4(nrow-1,2) = "-"
160 text4(nrow-1,3) = "-"
161 text4(nrow-1,4) = "-"
162 text4(nrow-1,5) = "-"
163 text4(nrow-1,6) = sprintf("%.2f",beta_4_avg)
164 text4(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
170 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
172 header = (/"<HTML>" \
174 ,"<TITLE>CLAMP metrics</TITLE>" \
181 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
183 ," <th bgcolor=DDDDDD >Station</th>" \
184 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
185 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
186 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
187 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
188 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
189 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
190 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
191 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
194 table_footer = "</table>"
198 lines = new(50000,string)
201 set_line(lines,nline,header)
202 set_line(lines,nline,table_header)
203 ;-----------------------------------------------
207 set_line(lines,nline,row_header)
219 set_line(lines,nline,"<th>"+txt1+"</th>")
220 set_line(lines,nline,"<th>"+txt2+"</th>")
221 set_line(lines,nline,"<th>"+txt3+"</th>")
222 set_line(lines,nline,"<th>"+txt4+"</th>")
223 set_line(lines,nline,"<th>"+txt5+"</th>")
224 set_line(lines,nline,"<th>"+txt6+"</th>")
225 set_line(lines,nline,"<th>"+txt7+"</th>")
226 set_line(lines,nline,"<th>"+txt8+"</th>")
227 set_line(lines,nline,"<th>"+txt9+"</th>")
229 set_line(lines,nline,row_footer)
231 ;-----------------------------------------------
232 set_line(lines,nline,table_footer)
233 set_line(lines,nline,footer)
235 ; Now write to an HTML file.
236 idx = ind(.not.ismissing(lines))
237 if(.not.any(ismissing(idx))) then
238 asciiwrite(output_html,lines(idx))
246 delete (table_header)
249 ;------------------------------------------------
250 ; read biome data: model
252 biome_name_mod = "Model PFT Class"
254 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
255 film = "class_pft_"+model_grid+".nc"
257 fm = addfile(dirm+film,"r")
259 classmod = fm->CLASS_PFT
263 ; model data has 17 land-type classes
267 ;------------------------------------------------
268 ; read biome data: observed
270 biome_name_ob = "MODIS LandCover"
272 diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
273 filo = "land_class_"+model_grid+".nc"
275 fo = addfile(diro+filo,"r")
277 classob = tofloat(fo->LAND_CLASS)
281 ; input observed data has 20 land-type classes
285 ;********************************************************************
286 ; use land-type class to bin the data in equally spaced ranges
287 ;********************************************************************
289 ; using observed biome class
291 ; using model biome class
295 range = fspan(0,nclassn-1,nclassn)
298 ; Use this range information to grab all the values in a
299 ; particular range, and then take an average.
303 xvalues = new((/2,nx/),float)
304 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
305 dx = xvalues(0,1) - xvalues(0,0) ; range width
306 dx4 = dx/4 ; 1/4 of the range
307 xvalues(1,:) = xvalues(0,:) - dx/5.
309 ;==============================
311 ;==============================
313 ; using observed biome class
314 ; base_1D = ndtooned(classob)
315 ; using model biome class
316 base_1D = ndtooned(classmod)
318 data1_1D = ndtooned(npp_i)
319 data2_1D = ndtooned(npp_f)
323 yvalues = new((/2,nx/),float)
324 count = new((/2,nx/),float)
328 ; See if we are doing data1 (nd=0) or data2 (nd=1).
338 ; Loop through each range, using base.
341 if (i.ne.(nr-2)) then
343 ; print("In range ["+range(i)+","+range(i+1)+")")
344 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
347 ; print("In range ["+range(i)+",)")
348 idx = ind(base.ge.range(i))
353 if(.not.any(ismissing(idx))) then
354 yvalues(nd,i) = avg(data(idx))
355 count(nd,i) = dimsizes(idx)
357 yvalues(nd,i) = yvalues@_FillValue
361 ;#############################################################
362 ;using observed biome class:
364 ; set the following 4 classes to _FillValue:
365 ; Water Bodies(0), Urban and Build-Up(13),
366 ; Permenant Snow and Ice(15), Unclassified(17)
368 ; if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
369 ; yvalues(nd,i) = yvalues@_FillValue
372 ;#############################################################
374 ;#############################################################
375 ;using model biome class:
377 ; set the following 4 classes to _FillValue:
378 ; (3)Needleleaf Deciduous Boreal Tree,
379 ; (8)Broadleaf Deciduous Boreal Tree,
380 ; (9)Broadleaf Evergreen Shrub,
383 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
384 yvalues(nd,i) = yvalues@_FillValue
387 ;#############################################################
389 ; print(nd + ": " + count(nd,i) + " points, avg = " + yvalues(nd,i))
391 ; Clean up for next time in loop.
399 ;============================
401 ;============================
408 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
412 uu_count = u_count(good)
413 vv_count = v_count(good)
415 n_biome = dimsizes(uu)
416 beta_biome = new((/n_biome/),float)
418 beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
420 ;beta_biome_avg = avg(beta_biome)
421 beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
423 ;print (beta_biome_avg)
425 ;===========================
426 ; for html table - biome
427 ;===========================
429 output_html = "table_biome.html"
431 ; column (not including header column)
433 col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
435 ncol = dimsizes(col_head)
437 ; row (not including header row)
439 ;----------------------------------------------------
440 ; using observed biome class:
442 ; row_head = (/"Evergreen Needleleaf Forests" \
443 ; ,"Evergreen Broadleaf Forests" \
444 ; ,"Deciduous Needleleaf Forest" \
445 ; ,"Deciduous Broadleaf Forests" \
447 ; ,"Closed Bushlands" \
448 ; ,"Open Bushlands" \
449 ; ,"Woody Savannas (S. Hem.)" \
450 ; ,"Savannas (S. Hem.)" \
452 ; ,"Permanent Wetlands" \
454 ; ,"Cropland/Natural Vegetation Mosaic" \
455 ; ,"Barren or Sparsely Vegetated" \
456 ; ,"Woody Savannas (N. Hem.)" \
457 ; ,"Savannas (N. Hem.)" \
461 ;----------------------------------------------------
462 ; using model biome class:
464 row_head = (/"Not Vegetated" \
465 ,"Needleleaf Evergreen Temperate Tree" \
466 ,"Needleleaf Evergreen Boreal Tree" \
467 ; ,"Needleleaf Deciduous Boreal Tree" \
468 ,"Broadleaf Evergreen Tropical Tree" \
469 ,"Broadleaf Evergreen Temperate Tree" \
470 ,"Broadleaf Deciduous Tropical Tree" \
471 ,"Broadleaf Deciduous Temperate Tree" \
472 ; ,"Broadleaf Deciduous Boreal Tree" \
473 ; ,"Broadleaf Evergreen Shrub" \
474 ,"Broadleaf Deciduous Temperate Shrub" \
475 ,"Broadleaf Deciduous Boreal Shrub" \
477 ,"C3 Non-Arctic Grass" \
484 nrow = dimsizes(row_head)
486 ; arrays to be passed to table.
487 text4 = new ((/nrow, ncol/),string )
490 text4(i,0) = sprintf("%.1f",co2_i)
491 text4(i,1) = sprintf("%.1f",co2_f)
492 text4(i,2) = sprintf("%.1f",uu(i))
493 text4(i,3) = sprintf("%.1f",vv(i))
494 text4(i,4) = sprintf("%.2f",beta_biome(i))
496 text4(nrow-1,0) = "-"
497 text4(nrow-1,1) = "-"
498 text4(nrow-1,2) = "-"
499 text4(nrow-1,3) = "-"
500 text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
502 ;**************************************************
504 ;**************************************************
506 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
508 header = (/"<HTML>" \
510 ,"<TITLE>CLAMP metrics</TITLE>" \
517 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
519 ," <th bgcolor=DDDDDD >Biome Class</th>" \
520 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
521 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
522 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
523 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
524 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
527 table_footer = "</table>"
531 lines = new(50000,string)
534 set_line(lines,nline,header)
535 set_line(lines,nline,table_header)
536 ;-----------------------------------------------
540 set_line(lines,nline,row_header)
549 set_line(lines,nline,"<th>"+txt1+"</th>")
550 set_line(lines,nline,"<th>"+txt2+"</th>")
551 set_line(lines,nline,"<th>"+txt3+"</th>")
552 set_line(lines,nline,"<th>"+txt4+"</th>")
553 set_line(lines,nline,"<th>"+txt5+"</th>")
554 set_line(lines,nline,"<th>"+txt6+"</th>")
556 set_line(lines,nline,row_footer)
558 ;-----------------------------------------------
559 set_line(lines,nline,table_footer)
560 set_line(lines,nline,footer)
562 ; Now write to an HTML file.
563 idx = ind(.not.ismissing(lines))
564 if(.not.any(ismissing(idx))) then
565 asciiwrite(output_html,lines(idx))