|
1 ;******************************************************** |
|
2 ; hardwired: co2_i = 283.1878 |
|
3 ; co2_f = 364.1252 |
|
4 ; |
|
5 ; beta_4_ob = 0.6 |
|
6 ;************************************************************** |
|
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
10 ;************************************************************** |
|
11 procedure set_line(lines:string,nline:integer,newlines:string) |
|
12 begin |
|
13 ; add line to ascci/html file |
|
14 |
|
15 nnewlines = dimsizes(newlines) |
|
16 if(nline+nnewlines-1.ge.dimsizes(lines)) |
|
17 print("set_line: bad index, not setting anything.") |
|
18 return |
|
19 end if |
|
20 lines(nline:nline+nnewlines-1) = newlines |
|
21 ; print ("lines = " + lines(nline:nline+nnewlines-1)) |
|
22 nline = nline + nnewlines |
|
23 return |
|
24 end |
|
25 ;************************************************************** |
|
26 ; Main code. |
|
27 begin |
|
28 |
|
29 plot_type = "ps" |
|
30 plot_type_new = "png" |
|
31 |
|
32 ;------------------------------------------------------ |
|
33 ; edit table.html of current model for movel1_vs_model2 |
|
34 |
|
35 if (isvar("compare")) then |
|
36 html_name2 = compare+"/table.html" |
|
37 html_new2 = html_name2 +".new" |
|
38 end if |
|
39 |
|
40 ;------------------------------------------------------- |
|
41 ; edit table.html for current model |
|
42 |
|
43 html_name = model_name+"/table.html" |
|
44 html_new = html_name +".new" |
|
45 |
|
46 ;------------------------------------------------------- |
|
47 ; read model data |
|
48 |
|
49 ;############################################################### |
|
50 ; hardwired for model data |
|
51 |
|
52 ; these values correspond to the start and the end of model data |
|
53 co2_i = 283.1878 |
|
54 co2_f = 364.1252 |
|
55 |
|
56 film_i = film6 |
|
57 film_f = film5 |
|
58 ;############################################################### |
|
59 |
|
60 fm_i = addfile (dirm+film_i,"r") |
|
61 fm_f = addfile (dirm+film_f,"r") |
|
62 |
|
63 xm = fm_f->lon |
|
64 ym = fm_f->lat |
|
65 |
|
66 npp_i = fm_i->NPP |
|
67 npp_f = fm_f->NPP |
|
68 |
|
69 delete (fm_i) |
|
70 delete (fm_f) |
|
71 |
|
72 ;Units for these variables are: |
|
73 ;npp_i: g C/m^2/s |
|
74 |
|
75 nsec_per_year = 60*60*24*365 |
|
76 |
|
77 npp_i = npp_i * nsec_per_year |
|
78 npp_f = npp_f * nsec_per_year |
|
79 |
|
80 unit = "gC/m2/year" |
|
81 |
|
82 ;------------------------------ |
|
83 ; get landfrac data |
|
84 |
|
85 film_l = "lnd_" + model_grid + ".nc" |
|
86 fm_l = addfile (dirs+film_l,"r") |
|
87 landfrac = fm_l->landfrac |
|
88 |
|
89 npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:) |
|
90 npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:) |
|
91 |
|
92 delete (fm_l) |
|
93 delete (landfrac) |
|
94 |
|
95 ;----------------------------- |
|
96 ; read biome data: model |
|
97 |
|
98 biome_name_mod = "Model PFT Class" |
|
99 |
|
100 film_c = "class_pft_"+model_grid+".nc" |
|
101 fm_c = addfile (dirs+film_c,"r") |
|
102 classmod = fm_c->CLASS_PFT |
|
103 |
|
104 delete (fm_c) |
|
105 |
|
106 ; model data has 17 land-type classes |
|
107 nclass_mod = 17 |
|
108 |
|
109 ;--------------------------------------------------- |
|
110 ; read data: observed at stations |
|
111 |
|
112 station = (/"DukeFACE" \ |
|
113 ,"AspenFACE" \ |
|
114 ,"ORNL-FACE" \ |
|
115 ,"POP-EUROFACE" \ |
|
116 /) |
|
117 |
|
118 lat_ob = (/ 35.58, 45.40, 35.54, 42.22/) |
|
119 lon_ob = (/-79.05, -89.37, -84.20, 11.48/) |
|
120 lon_obx = where(lon_ob.lt.0.,lon_ob+360.,lon_ob) |
|
121 |
|
122 n_sta = dimsizes(station) |
|
123 beta_4_ob = new((/n_sta/),float) |
|
124 |
|
125 ;################################################### |
|
126 ; this is a hardwired value |
|
127 beta_4_ob = 0.60 |
|
128 ;################################################### |
|
129 ;--------------------------------------------------- |
|
130 ; get model data at station |
|
131 |
|
132 npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_obx,lat_ob,0) |
|
133 |
|
134 npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_obx,lat_ob,0) |
|
135 |
|
136 ;--------------------------------------------------- |
|
137 ;compute beta_4 |
|
138 |
|
139 score_max = 3. |
|
140 |
|
141 beta_4 = new((/n_sta/),float) |
|
142 |
|
143 beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i) |
|
144 |
|
145 beta_4_avg = avg(beta_4) |
|
146 |
|
147 bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob))) |
|
148 Mbeta = (1. - (bias/n_sta))*score_max |
|
149 M_beta = sprintf("%.2f", Mbeta) |
|
150 |
|
151 ;========================= |
|
152 ; for html table - station |
|
153 ;========================= |
|
154 |
|
155 output_html = "table_station.html" |
|
156 |
|
157 ; column (not including header column) |
|
158 |
|
159 col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/) |
|
160 |
|
161 ncol = dimsizes(col_head) |
|
162 |
|
163 ; row (not including header row) |
|
164 row_head = (/"DukeFACE" \ |
|
165 ,"AspenFACE" \ |
|
166 ,"ORNL-FACE" \ |
|
167 ,"POP-EUROFACE" \ |
|
168 ,"All Station" \ |
|
169 /) |
|
170 nrow = dimsizes(row_head) |
|
171 |
|
172 ; arrays to be passed to table. |
|
173 text = new ((/nrow, ncol/),string ) |
|
174 |
|
175 do i=0,nrow-2 |
|
176 text(i,0) = sprintf("%.1f",lat_ob(i)) |
|
177 text(i,1) = sprintf("%.1f",lon_ob(i)) |
|
178 text(i,2) = sprintf("%.1f",co2_i) |
|
179 text(i,3) = sprintf("%.1f",co2_f) |
|
180 text(i,4) = sprintf("%.1f",npp_i_4(0,i)) |
|
181 text(i,5) = sprintf("%.1f",npp_f_4(0,i)) |
|
182 text(i,6) = sprintf("%.2f",beta_4(i)) |
|
183 text(i,7) = "-" |
|
184 end do |
|
185 text(nrow-1,0) = "-" |
|
186 text(nrow-1,1) = "-" |
|
187 text(nrow-1,2) = "-" |
|
188 text(nrow-1,3) = "-" |
|
189 text(nrow-1,4) = "-" |
|
190 text(nrow-1,5) = "-" |
|
191 text(nrow-1,6) = sprintf("%.2f",beta_4_avg) |
|
192 text(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob)) |
|
193 |
|
194 ;----------- |
|
195 ; html table |
|
196 ;----------- |
|
197 |
|
198 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>" |
|
199 |
|
200 header = (/"<HTML>" \ |
|
201 ,"<HEAD>" \ |
|
202 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
203 ,"</HEAD>" \ |
|
204 ,header_text \ |
|
205 /) |
|
206 footer = "</HTML>" |
|
207 |
|
208 table_header = (/ \ |
|
209 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \ |
|
210 ,"<tr>" \ |
|
211 ," <th bgcolor=DDDDDD >Station</th>" \ |
|
212 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \ |
|
213 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \ |
|
214 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \ |
|
215 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \ |
|
216 ," <th bgcolor=DDDDDD >"+col_head(4)+"<br>("+unit+")</th>" \ |
|
217 ," <th bgcolor=DDDDDD >"+col_head(5)+"<br>("+unit+")</th>" \ |
|
218 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \ |
|
219 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \ |
|
220 ,"</tr>" \ |
|
221 /) |
|
222 table_footer = "</table>" |
|
223 row_header = "<tr>" |
|
224 row_footer = "</tr>" |
|
225 |
|
226 lines = new(50000,string) |
|
227 nline = 0 |
|
228 |
|
229 set_line(lines,nline,header) |
|
230 set_line(lines,nline,table_header) |
|
231 ;----------------------------------------------- |
|
232 ;row of table |
|
233 |
|
234 do n = 0,nrow-1 |
|
235 set_line(lines,nline,row_header) |
|
236 |
|
237 txt1 = row_head(n) |
|
238 txt2 = text(n,0) |
|
239 txt3 = text(n,1) |
|
240 txt4 = text(n,2) |
|
241 txt5 = text(n,3) |
|
242 txt6 = text(n,4) |
|
243 txt7 = text(n,5) |
|
244 txt8 = text(n,6) |
|
245 txt9 = text(n,7) |
|
246 |
|
247 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
248 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
249 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
250 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
251 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
252 set_line(lines,nline,"<th>"+txt6+"</th>") |
|
253 set_line(lines,nline,"<th>"+txt7+"</th>") |
|
254 set_line(lines,nline,"<th>"+txt8+"</th>") |
|
255 set_line(lines,nline,"<th>"+txt9+"</th>") |
|
256 |
|
257 set_line(lines,nline,row_footer) |
|
258 end do |
|
259 ;----------------------------------------------- |
|
260 set_line(lines,nline,table_footer) |
|
261 set_line(lines,nline,footer) |
|
262 |
|
263 ; Now write to an HTML file. |
|
264 idx = ind(.not.ismissing(lines)) |
|
265 if(.not.any(ismissing(idx))) then |
|
266 asciiwrite(output_html,lines(idx)) |
|
267 else |
|
268 print ("error?") |
|
269 end if |
|
270 |
|
271 delete (col_head) |
|
272 delete (row_head) |
|
273 delete (text) |
|
274 delete (table_header) |
|
275 delete (idx) |
|
276 |
|
277 ;******************************************************************** |
|
278 ; use land-type class to bin the data in equally spaced ranges |
|
279 ;******************************************************************** |
|
280 |
|
281 ; using model biome class |
|
282 nclass = nclass_mod |
|
283 |
|
284 range = fspan(0,nclass,nclass+1) |
|
285 |
|
286 ; Use this range information to grab all the values in a |
|
287 ; particular range, and then take an average. |
|
288 |
|
289 nx = dimsizes(range) - 1 |
|
290 |
|
291 ;============================== |
|
292 ; put data into bins |
|
293 ;============================== |
|
294 |
|
295 ; for model data and observed |
|
296 data_n = 2 |
|
297 |
|
298 ; using model biome class |
|
299 |
|
300 base = ndtooned(classmod) |
|
301 |
|
302 ; output |
|
303 |
|
304 yvalues = new((/data_n,nx/),float) |
|
305 count = new((/data_n,nx/),float) |
|
306 |
|
307 ; Loop through each range, using base |
|
308 |
|
309 do i=0,nx-1 |
|
310 |
|
311 if (i.ne.(nx-1)) then |
|
312 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) |
|
313 else |
|
314 idx = ind(base.ge.range(i)) |
|
315 end if |
|
316 |
|
317 ; loop through each dataset |
|
318 |
|
319 do n = 0,data_n-1 |
|
320 |
|
321 if (n .eq. 0) then |
|
322 data = ndtooned(npp_i) |
|
323 end if |
|
324 |
|
325 if (n .eq. 1) then |
|
326 data = ndtooned(npp_f) |
|
327 end if |
|
328 |
|
329 ; Calculate average |
|
330 |
|
331 if (.not.any(ismissing(idx))) then |
|
332 yvalues(n,i) = avg(data(idx)) |
|
333 count(n,i) = dimsizes(idx) |
|
334 else |
|
335 yvalues(n,i) = yvalues@_FillValue |
|
336 count(n,i) = 0 |
|
337 end if |
|
338 |
|
339 ;############################################################# |
|
340 ; using model biome class: |
|
341 ; |
|
342 ; set the following 4 classes to _FillValue: |
|
343 ; (3)Needleleaf Deciduous Boreal Tree, |
|
344 ; (8)Broadleaf Deciduous Boreal Tree, |
|
345 ; (9)Broadleaf Evergreen Shrub, |
|
346 ; (16)Wheat |
|
347 |
|
348 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then |
|
349 yvalues(n,i) = yvalues@_FillValue |
|
350 count(n,i) = 0 |
|
351 end if |
|
352 ;############################################################# |
|
353 |
|
354 delete(data) |
|
355 end do ; n-loop |
|
356 |
|
357 delete(idx) |
|
358 end do ; i-loop |
|
359 |
|
360 delete (base) |
|
361 delete (npp_i) |
|
362 delete (npp_f) |
|
363 |
|
364 ;============================ |
|
365 ;compute beta |
|
366 ;============================ |
|
367 |
|
368 u = yvalues(0,:) |
|
369 v = yvalues(1,:) |
|
370 u_count = count(0,:) |
|
371 v_count = count(1,:) |
|
372 |
|
373 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
374 |
|
375 uu = u(good) |
|
376 vv = v(good) |
|
377 uu_count = u_count(good) |
|
378 vv_count = v_count(good) |
|
379 |
|
380 n_biome = dimsizes(uu) |
|
381 beta_biome = new((/n_biome/),float) |
|
382 |
|
383 beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i) |
|
384 |
|
385 beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i) |
|
386 |
|
387 ;=========================== |
|
388 ; for html table - biome |
|
389 ;=========================== |
|
390 |
|
391 output_html = "table_biome.html" |
|
392 |
|
393 ; column (not including header column) |
|
394 |
|
395 col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/) |
|
396 |
|
397 ncol = dimsizes(col_head) |
|
398 |
|
399 ; row (not including header row) |
|
400 |
|
401 ;---------------------------------------------------- |
|
402 ; using model biome class: |
|
403 ; |
|
404 row_head = (/"Not Vegetated" \ |
|
405 ,"Needleleaf Evergreen Temperate Tree" \ |
|
406 ,"Needleleaf Evergreen Boreal Tree" \ |
|
407 ; ,"Needleleaf Deciduous Boreal Tree" \ |
|
408 ,"Broadleaf Evergreen Tropical Tree" \ |
|
409 ,"Broadleaf Evergreen Temperate Tree" \ |
|
410 ,"Broadleaf Deciduous Tropical Tree" \ |
|
411 ,"Broadleaf Deciduous Temperate Tree" \ |
|
412 ; ,"Broadleaf Deciduous Boreal Tree" \ |
|
413 ; ,"Broadleaf Evergreen Shrub" \ |
|
414 ,"Broadleaf Deciduous Temperate Shrub" \ |
|
415 ,"Broadleaf Deciduous Boreal Shrub" \ |
|
416 ,"C3 Arctic Grass" \ |
|
417 ,"C3 Non-Arctic Grass" \ |
|
418 ,"C4 Grass" \ |
|
419 ,"Corn" \ |
|
420 ; ,"Wheat" \ |
|
421 ,"All Biome" \ |
|
422 /) |
|
423 |
|
424 nrow = dimsizes(row_head) |
|
425 |
|
426 ; arrays to be passed to table. |
|
427 text = new ((/nrow, ncol/),string ) |
|
428 |
|
429 do i=0,nrow-2 |
|
430 text(i,0) = sprintf("%.1f",co2_i) |
|
431 text(i,1) = sprintf("%.1f",co2_f) |
|
432 text(i,2) = sprintf("%.1f",uu(i)) |
|
433 text(i,3) = sprintf("%.1f",vv(i)) |
|
434 text(i,4) = sprintf("%.2f",beta_biome(i)) |
|
435 end do |
|
436 text(nrow-1,0) = "-" |
|
437 text(nrow-1,1) = "-" |
|
438 text(nrow-1,2) = "-" |
|
439 text(nrow-1,3) = "-" |
|
440 text(nrow-1,4) = sprintf("%.2f",beta_biome_avg) |
|
441 |
|
442 ;************************************************** |
|
443 ; html table |
|
444 ;************************************************** |
|
445 |
|
446 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>" |
|
447 |
|
448 header = (/"<HTML>" \ |
|
449 ,"<HEAD>" \ |
|
450 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
451 ,"</HEAD>" \ |
|
452 ,header_text \ |
|
453 /) |
|
454 footer = "</HTML>" |
|
455 |
|
456 table_header = (/ \ |
|
457 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \ |
|
458 ,"<tr>" \ |
|
459 ," <th bgcolor=DDDDDD >Biome Class</th>" \ |
|
460 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \ |
|
461 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \ |
|
462 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \ |
|
463 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \ |
|
464 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \ |
|
465 ,"</tr>" \ |
|
466 /) |
|
467 table_footer = "</table>" |
|
468 row_header = "<tr>" |
|
469 row_footer = "</tr>" |
|
470 |
|
471 lines = new(50000,string) |
|
472 nline = 0 |
|
473 |
|
474 set_line(lines,nline,header) |
|
475 set_line(lines,nline,table_header) |
|
476 ;----------------------------------------------- |
|
477 ;row of table |
|
478 |
|
479 do n = 0,nrow-1 |
|
480 set_line(lines,nline,row_header) |
|
481 |
|
482 txt1 = row_head(n) |
|
483 txt2 = text(n,0) |
|
484 txt3 = text(n,1) |
|
485 txt4 = text(n,2) |
|
486 txt5 = text(n,3) |
|
487 txt6 = text(n,4) |
|
488 |
|
489 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
490 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
491 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
492 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
493 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
494 set_line(lines,nline,"<th>"+txt6+"</th>") |
|
495 |
|
496 set_line(lines,nline,row_footer) |
|
497 end do |
|
498 ;----------------------------------------------- |
|
499 set_line(lines,nline,table_footer) |
|
500 set_line(lines,nline,footer) |
|
501 |
|
502 ; Now write to an HTML file |
|
503 |
|
504 idx = ind(.not.ismissing(lines)) |
|
505 if(.not.any(ismissing(idx))) then |
|
506 asciiwrite(output_html,lines(idx)) |
|
507 else |
|
508 print ("error?") |
|
509 end if |
|
510 delete (idx) |
|
511 |
|
512 ;************************************************************************************** |
|
513 ; update score |
|
514 ;************************************************************************************** |
|
515 if (isvar("compare")) then |
|
516 system("sed -e '1,/M_beta/s/M_beta/"+M_beta+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
517 "mv -f "+html_new2+" "+html_name2) |
|
518 end if |
|
519 |
|
520 system("sed s#M_beta#"+M_beta+"# "+html_name+" > "+html_new+";"+ \ |
|
521 "mv -f "+html_new+" "+html_name) |
|
522 |
|
523 ;*************************************************************************** |
|
524 ; get total score and write to file |
|
525 ;*************************************************************************** |
|
526 M_total = Mbeta |
|
527 |
|
528 asciiwrite("M_save.beta", M_total) |
|
529 |
|
530 delete (M_total) |
|
531 |
|
532 ;*************************************************************************** |
|
533 ; output plot and html |
|
534 ;*************************************************************************** |
|
535 output_dir = model_name+"/beta" |
|
536 |
|
537 system("mv *.html " + output_dir) |
|
538 ;*************************************************************************** |
|
539 |
|
540 end |
|
541 |