|
1 ;************************************************************** |
|
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
5 ;************************************************************** |
|
6 procedure set_line(lines:string,nline:integer,newlines:string) |
|
7 begin |
|
8 ; add line to ascci/html file |
|
9 |
|
10 nnewlines = dimsizes(newlines) |
|
11 if(nline+nnewlines-1.ge.dimsizes(lines)) |
|
12 print("set_line: bad index, not setting anything.") |
|
13 return |
|
14 end if |
|
15 lines(nline:nline+nnewlines-1) = newlines |
|
16 ; print ("lines = " + lines(nline:nline+nnewlines-1)) |
|
17 nline = nline + nnewlines |
|
18 return |
|
19 end |
|
20 ;************************************************************** |
|
21 ; Main code. |
|
22 begin |
|
23 |
|
24 plot_type = "ps" |
|
25 plot_type_new = "png" |
|
26 |
|
27 ;---------------------------------------------------------- |
|
28 ; edit current model for movel1_vs_model2 |
|
29 |
|
30 if (isvar("compare")) then |
|
31 html_name2 = compare+"/table.html" |
|
32 html_new2 = html_name2 +".new" |
|
33 end if |
|
34 |
|
35 ;---------------------------------------------------------- |
|
36 ; edit table.html for current model |
|
37 |
|
38 html_name = model_name+"/table.html" |
|
39 html_new = html_name +".new" |
|
40 |
|
41 ;---------------------------------------------------------- |
|
42 ; get biome data: model |
|
43 |
|
44 biome_name_mod = "Model PFT Class" |
|
45 |
|
46 film_c = "class_pft_"+ model_grid +".nc" |
|
47 fm_c = addfile (dirs+film_c,"r") |
|
48 classmod = fm_c->CLASS_PFT |
|
49 |
|
50 delete (fm_c) |
|
51 |
|
52 ; model data has 17 land-type classes |
|
53 nclass_mod = 17 |
|
54 |
|
55 ;------------------------------------ |
|
56 ; get model landfrac and area |
|
57 |
|
58 film_l = "lnd_"+ model_grid +".nc" |
|
59 fm_l = addfile (dirs+film_l,"r") |
|
60 landfrac = fm_l->landfrac |
|
61 area = fm_l->area |
|
62 |
|
63 delete (fm_l) |
|
64 |
|
65 ; change area from km**2 to m**2 |
|
66 area = area * 1.e6 |
|
67 |
|
68 ;------------------------------------- |
|
69 ; take into account landfrac |
|
70 |
|
71 ; area = area * landfrac |
|
72 |
|
73 ; delete (landfrac) |
|
74 |
|
75 ;--------------------------------------------------- |
|
76 ; read data: model, group 1 |
|
77 |
|
78 fm = addfile (dirm+film4,"r") |
|
79 |
|
80 NPP1 = fm->NPP |
|
81 |
|
82 leafc = fm->LEAFC |
|
83 woodc = fm->WOODC |
|
84 frootc = fm->FROOTC |
|
85 VegC = leafc |
|
86 VegC = leafc + woodc + frootc |
|
87 |
|
88 litterc = fm->LITTERC |
|
89 cwdc = fm->CWDC |
|
90 LiCwC = litterc |
|
91 LiCwC = litterc + cwdc |
|
92 |
|
93 SoilC = fm->SOILC |
|
94 |
|
95 delete (fm) |
|
96 ;--------------------------------------------------- |
|
97 ; read data: model, group 2 |
|
98 |
|
99 fm = addfile (dirm+film5,"r") |
|
100 |
|
101 NPP2 = fm->NPP |
|
102 NEE2 = fm->NEE |
|
103 GPP2 = fm->GPP |
|
104 |
|
105 delete (fm) |
|
106 ;--------------------------------------------------- |
|
107 ; Units for these variables are: |
|
108 |
|
109 ;NPP1: g C/m^2/s |
|
110 ;NPP2: g C/m^2/s |
|
111 ;NEE2: g C/m^2/s |
|
112 ;GPP2: g C/m^2/s |
|
113 |
|
114 ;VegC: g C/m^2 |
|
115 ;LiCwC: g C/m^2 |
|
116 ;SoilC: g C/m^2 |
|
117 |
|
118 nsec_per_year = 60*60*24*365 |
|
119 |
|
120 ; change unit to g C/m^2/year |
|
121 |
|
122 NPP1 = NPP1 * nsec_per_year * conform(NPP1,landfrac,(/1,2/)) |
|
123 NPP2 = NPP2 * nsec_per_year * conform(NPP2,landfrac,(/1,2/)) |
|
124 NEE2 = NEE2 * nsec_per_year * conform(NEE2,landfrac,(/1,2/)) |
|
125 GPP2 = GPP2 * nsec_per_year * conform(GPP2,landfrac,(/1,2/)) |
|
126 |
|
127 VegC = VegC * conform(VegC,landfrac,(/1,2/)) |
|
128 LiCwC = LiCwC * conform(LiCwC,landfrac,(/1,2/)) |
|
129 SoilC = SoilC * conform(SoilC,landfrac,(/1,2/)) |
|
130 |
|
131 data_n = 8 |
|
132 |
|
133 ;******************************************************************* |
|
134 ; Calculate "nice" bins for binning the data in equally spaced ranges |
|
135 ;******************************************************************** |
|
136 |
|
137 ; using model biome class |
|
138 nclass = nclass_mod |
|
139 |
|
140 range = fspan(0,nclass,nclass+1) |
|
141 |
|
142 ; print (range) |
|
143 ; Use this range information to grab all the values in a |
|
144 ; particular range, and then take an average. |
|
145 |
|
146 nx = dimsizes(range) - 1 |
|
147 |
|
148 ;============================== |
|
149 ; put data into bins |
|
150 ;============================== |
|
151 |
|
152 ; using observed biome class |
|
153 ; base = ndtooned(classob) |
|
154 ; using model biome class |
|
155 base = ndtooned(classmod) |
|
156 |
|
157 area_1d = ndtooned(area) |
|
158 |
|
159 ; output |
|
160 |
|
161 yvalues = new((/data_n,nx/),float) ; (per m2) |
|
162 yvalues_t = new((/data_n,nx/),float) ; (per biome) |
|
163 |
|
164 ; Loop through each range, using base. |
|
165 |
|
166 do i=0,nx-1 |
|
167 |
|
168 if (i.ne.(nx-1)) then |
|
169 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) |
|
170 else |
|
171 idx = ind(base.ge.range(i)) |
|
172 end if |
|
173 |
|
174 do n = 0,data_n-1 |
|
175 |
|
176 if (n.eq.0) then |
|
177 data = ndtooned(area) |
|
178 end if |
|
179 |
|
180 if (n.eq.1) then |
|
181 data = ndtooned(NPP1) |
|
182 end if |
|
183 |
|
184 if (n.eq.2) then |
|
185 data = ndtooned(VegC) |
|
186 end if |
|
187 |
|
188 if (n.eq.3) then |
|
189 data = ndtooned(LiCwC) |
|
190 end if |
|
191 |
|
192 if (n.eq.4) then |
|
193 data = ndtooned(SoilC) |
|
194 end if |
|
195 |
|
196 if (n.eq.5) then |
|
197 data = ndtooned(NPP2) |
|
198 end if |
|
199 |
|
200 if (n.eq.6) then |
|
201 data = ndtooned(NEE2) |
|
202 end if |
|
203 |
|
204 if (n.eq.7) then |
|
205 data = ndtooned(GPP2) |
|
206 end if |
|
207 |
|
208 ; Calculate sum and average |
|
209 |
|
210 if (.not.any(ismissing(idx))) then |
|
211 if (n.eq.0) then |
|
212 yvalues(n,i) = sum(data(idx)) |
|
213 yvalues_t(n,i) = sum(data(idx)) |
|
214 else |
|
215 yvalues(n,i) = avg(data(idx)) |
|
216 yvalues_t(n,i) = sum(data(idx)*area_1d(idx)) |
|
217 end if |
|
218 else |
|
219 yvalues(n,i) = yvalues@_FillValue |
|
220 yvalues_t(n,i) = yvalues@_FillValue |
|
221 end if |
|
222 |
|
223 ;############################################################# |
|
224 ; using model biome class: |
|
225 ; set the following 4 classes to _FillValue: |
|
226 ; (3)Needleleaf Deciduous Boreal Tree, |
|
227 ; (8)Broadleaf Deciduous Boreal Tree, |
|
228 ; (9)Broadleaf Evergreen Shrub, |
|
229 ; (16)Wheat |
|
230 |
|
231 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then |
|
232 yvalues(n,i) = yvalues@_FillValue |
|
233 yvalues_t(n,i) = yvalues@_FillValue |
|
234 end if |
|
235 ;############################################################# |
|
236 |
|
237 delete (data) |
|
238 end do |
|
239 |
|
240 delete (idx) |
|
241 end do |
|
242 |
|
243 delete (base) |
|
244 delete (area) |
|
245 delete (NPP1) |
|
246 delete (VegC) |
|
247 delete (LiCwC) |
|
248 delete (SoilC) |
|
249 delete (NPP2) |
|
250 delete (NEE2) |
|
251 delete (GPP2) |
|
252 |
|
253 ;---------------------------------------------------------------- |
|
254 ; data for table1 |
|
255 |
|
256 good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:))) |
|
257 ;print (good) |
|
258 |
|
259 area_g = yvalues(0,good) |
|
260 NPP1_g = yvalues(1,good) |
|
261 VegC_g = yvalues(2,good) |
|
262 LiCwC_g = yvalues(3,good) |
|
263 SoilC_g = yvalues(4,good) |
|
264 NPP2_g = yvalues(5,good) |
|
265 NEE2_g = yvalues(6,good) |
|
266 GPP2_g = yvalues(7,good) |
|
267 |
|
268 NPP_ratio = NPP2_g/NPP1_g |
|
269 |
|
270 n_biome = dimsizes(NPP1_g) |
|
271 |
|
272 ;----------------------------------------------------------------- |
|
273 ; data for table2 |
|
274 |
|
275 ; change unit from g to Pg (Peta gram) |
|
276 factor_unit = 1.e-15 |
|
277 |
|
278 NPP1_t = yvalues_t(1,good) * factor_unit |
|
279 VegC_t = yvalues_t(2,good) * factor_unit |
|
280 LiCwC_t = yvalues_t(3,good) * factor_unit |
|
281 SoilC_t = yvalues_t(4,good) * factor_unit |
|
282 NEE2_t = yvalues_t(6,good) * factor_unit |
|
283 GPP2_t = yvalues_t(7,good) * factor_unit |
|
284 |
|
285 delete (yvalues) |
|
286 delete (yvalues_t) |
|
287 |
|
288 ;------------------------------------------------------------- |
|
289 ; html table1 data |
|
290 |
|
291 ; column (not including header column) |
|
292 |
|
293 col_head = (/"Area (1.e12m2)" \ |
|
294 ,"NPP (gC/m2/yr)" \ |
|
295 ,"VegC (gC/m2)" \ |
|
296 ,"Litter+CWD (gC/m2)" \ |
|
297 ,"SoilC (gC/m2)" \ |
|
298 ,"NPP_ratio" \ |
|
299 ,"NEE (gC/m2/yr)" \ |
|
300 ,"GPP (gC/m2/yr)" \ |
|
301 /) |
|
302 |
|
303 ncol = dimsizes(col_head) |
|
304 |
|
305 ; row (not including header row) |
|
306 |
|
307 ; using model biome class: |
|
308 row_head = (/"Not Vegetated" \ |
|
309 ,"Needleleaf Evergreen Temperate Tree" \ |
|
310 ,"Needleleaf Evergreen Boreal Tree" \ |
|
311 ; ,"Needleleaf Deciduous Boreal Tree" \ |
|
312 ,"Broadleaf Evergreen Tropical Tree" \ |
|
313 ,"Broadleaf Evergreen Temperate Tree" \ |
|
314 ,"Broadleaf Deciduous Tropical Tree" \ |
|
315 ,"Broadleaf Deciduous Temperate Tree" \ |
|
316 ; ,"Broadleaf Deciduous Boreal Tree" \ |
|
317 ; ,"Broadleaf Evergreen Shrub" \ |
|
318 ,"Broadleaf Deciduous Temperate Shrub" \ |
|
319 ,"Broadleaf Deciduous Boreal Shrub" \ |
|
320 ,"C3 Arctic Grass" \ |
|
321 ,"C3 Non-Arctic Grass" \ |
|
322 ,"C4 Grass" \ |
|
323 ,"Corn" \ |
|
324 ; ,"Wheat" \ |
|
325 ,"All Biome" \ |
|
326 /) |
|
327 nrow = dimsizes(row_head) |
|
328 |
|
329 ; arrays to be passed to table. |
|
330 text = new ((/nrow, ncol/),string ) |
|
331 |
|
332 do i=0,nrow-2 |
|
333 text(i,0) = sprintf("%.1f",area_g(i)*1.e-12) |
|
334 text(i,1) = sprintf("%.1f",NPP1_g(i)) |
|
335 text(i,2) = sprintf("%.1f",VegC_g(i)) |
|
336 text(i,3) = sprintf("%.1f",LiCwC_g(i)) |
|
337 text(i,4) = sprintf("%.1f",SoilC_g(i)) |
|
338 text(i,5) = sprintf("%.2f",NPP_ratio(i)) |
|
339 text(i,6) = sprintf("%.1f",NEE2_g(i)) |
|
340 text(i,7) = sprintf("%.1f",GPP2_g(i)) |
|
341 end do |
|
342 |
|
343 ;------------------------------------------------------- |
|
344 ; create html table1 |
|
345 |
|
346 header_text = "<H1>NEE and Carbon Stocks and Fluxes: Model "+model_name+"</H1>" |
|
347 |
|
348 header = (/"<HTML>" \ |
|
349 ,"<HEAD>" \ |
|
350 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
351 ,"</HEAD>" \ |
|
352 ,header_text \ |
|
353 /) |
|
354 footer = "</HTML>" |
|
355 |
|
356 table_header = (/ \ |
|
357 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
|
358 ,"<tr>" \ |
|
359 ," <th bgcolor=DDDDDD >Biome Type</th>" \ |
|
360 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \ |
|
361 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \ |
|
362 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \ |
|
363 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \ |
|
364 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \ |
|
365 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \ |
|
366 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \ |
|
367 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \ |
|
368 ,"</tr>" \ |
|
369 /) |
|
370 table_footer = "</table>" |
|
371 row_header = "<tr>" |
|
372 row_footer = "</tr>" |
|
373 |
|
374 lines = new(50000,string) |
|
375 nline = 0 |
|
376 |
|
377 set_line(lines,nline,header) |
|
378 set_line(lines,nline,table_header) |
|
379 |
|
380 ;---------------------------- |
|
381 ;row of table |
|
382 |
|
383 do n = 0,nrow-2 |
|
384 set_line(lines,nline,row_header) |
|
385 |
|
386 txt0 = row_head(n) |
|
387 txt1 = text(n,0) |
|
388 txt2 = text(n,1) |
|
389 txt3 = text(n,2) |
|
390 txt4 = text(n,3) |
|
391 txt5 = text(n,4) |
|
392 txt6 = text(n,5) |
|
393 txt7 = text(n,6) |
|
394 txt8 = text(n,7) |
|
395 |
|
396 set_line(lines,nline,"<th>"+txt0+"</th>") |
|
397 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
398 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
399 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
400 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
401 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
402 set_line(lines,nline,"<th>"+txt6+"</th>") |
|
403 set_line(lines,nline,"<th>"+txt7+"</th>") |
|
404 set_line(lines,nline,"<th>"+txt8+"</th>") |
|
405 |
|
406 set_line(lines,nline,row_footer) |
|
407 end do |
|
408 ;---------------------------- |
|
409 set_line(lines,nline,table_footer) |
|
410 set_line(lines,nline,footer) |
|
411 |
|
412 ; Now write to an HTML file. |
|
413 |
|
414 output_html = "table_per_m2.html" |
|
415 |
|
416 idx = ind(.not.ismissing(lines)) |
|
417 if(.not.any(ismissing(idx))) then |
|
418 asciiwrite(output_html,lines(idx)) |
|
419 else |
|
420 print ("error?") |
|
421 end if |
|
422 |
|
423 delete (idx) |
|
424 |
|
425 delete (col_head) |
|
426 delete (row_head) |
|
427 delete (text) |
|
428 delete (table_header) |
|
429 |
|
430 ;----------------------------------------------------------------- |
|
431 ; html table2 data |
|
432 |
|
433 ; column (not including header column) |
|
434 |
|
435 col_head = (/"NPP (PgC/yr)" \ |
|
436 ,"VegC (PgC)" \ |
|
437 ,"Litter+CWD (PgC)" \ |
|
438 ,"SoilC (PgC)" \ |
|
439 ,"NEE (PgC/yr)" \ |
|
440 ,"GPP (PgC/yr)" \ |
|
441 ,"NPP timeseries" \ |
|
442 ,"NEE timeseries" \ |
|
443 ,"Fire timeseries" \ |
|
444 /) |
|
445 |
|
446 ncol = dimsizes(col_head) |
|
447 |
|
448 ; row (not including header row) |
|
449 |
|
450 ; using model biome class: |
|
451 row_head = (/"Not Vegetated" \ |
|
452 ,"Needleleaf Evergreen Temperate Tree" \ |
|
453 ,"Needleleaf Evergreen Boreal Tree" \ |
|
454 ; ,"Needleleaf Deciduous Boreal Tree" \ |
|
455 ,"Broadleaf Evergreen Tropical Tree" \ |
|
456 ,"Broadleaf Evergreen Temperate Tree" \ |
|
457 ,"Broadleaf Deciduous Tropical Tree" \ |
|
458 ,"Broadleaf Deciduous Temperate Tree" \ |
|
459 ; ,"Broadleaf Deciduous Boreal Tree" \ |
|
460 ; ,"Broadleaf Evergreen Shrub" \ |
|
461 ,"Broadleaf Deciduous Temperate Shrub" \ |
|
462 ,"Broadleaf Deciduous Boreal Shrub" \ |
|
463 ,"C3 Arctic Grass" \ |
|
464 ,"C3 Non-Arctic Grass" \ |
|
465 ,"C4 Grass" \ |
|
466 ,"Corn" \ |
|
467 ; ,"Wheat" \ |
|
468 ,"All Biome" \ |
|
469 /) |
|
470 nrow = dimsizes(row_head) |
|
471 |
|
472 ; arrays to be passed to table. |
|
473 text = new ((/nrow, ncol/),string ) |
|
474 |
|
475 do i=0,nrow-2 |
|
476 text(i,0) = sprintf("%.1f",NPP1_t(i)) |
|
477 text(i,1) = sprintf("%.1f",VegC_t(i)) |
|
478 text(i,2) = sprintf("%.1f",LiCwC_t(i)) |
|
479 text(i,3) = sprintf("%.1f",SoilC_t(i)) |
|
480 text(i,4) = sprintf("%.1f",NEE2_t(i)) |
|
481 text(i,5) = sprintf("%.1f",GPP2_t(i)) |
|
482 text(i,6) = "<a href=./NPP_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NPP_annual_biome_"+i+".png>annual_plot</a>" |
|
483 text(i,7) = "<a href=./NEE_monthly_biome_"+i+".png>monthly_plot</a> <br> <a href=./NEE_annual_biome_"+i+".png>annual_plot</a>" |
|
484 text(i,8) = "--" |
|
485 end do |
|
486 |
|
487 text(nrow-1,0) = sprintf("%.1f",sum(NPP1_t)) |
|
488 text(nrow-1,1) = sprintf("%.1f",sum(VegC_t)) |
|
489 text(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t)) |
|
490 text(nrow-1,3) = sprintf("%.1f",sum(SoilC_t)) |
|
491 text(nrow-1,4) = sprintf("%.1f",sum(NEE2_t)) |
|
492 text(nrow-1,5) = sprintf("%.1f",sum(GPP2_t)) |
|
493 text(nrow-1,6) = "<a href=./NPP_monthly_global.png>monthly_plot</a> <br> <a href=./NPP_annual_global.png>annual_plot</a>" |
|
494 text(nrow-1,7) = "<a href=./NEE_monthly_global.png>monthly_plot</a> <br> <a href=./NEE_annual_global.png>annual_plot</a>" |
|
495 text(nrow-1,8) = "--" |
|
496 |
|
497 ;************************************************** |
|
498 ; create html table2 |
|
499 ;************************************************** |
|
500 |
|
501 header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome): Model "+model_name+"</H1>" |
|
502 |
|
503 header = (/"<HTML>" \ |
|
504 ,"<HEAD>" \ |
|
505 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
506 ,"</HEAD>" \ |
|
507 ,header_text \ |
|
508 /) |
|
509 footer = "</HTML>" |
|
510 |
|
511 table_header = (/ \ |
|
512 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
|
513 ,"<tr>" \ |
|
514 ," <th bgcolor=DDDDDD >Biome Type</th>" \ |
|
515 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \ |
|
516 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \ |
|
517 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \ |
|
518 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \ |
|
519 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \ |
|
520 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \ |
|
521 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \ |
|
522 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \ |
|
523 ," <th bgcolor=DDDDDD >"+col_head(8)+"</th>" \ |
|
524 ,"</tr>" \ |
|
525 /) |
|
526 table_footer = "</table>" |
|
527 row_header = "<tr>" |
|
528 row_footer = "</tr>" |
|
529 |
|
530 lines = new(50000,string) |
|
531 nline = 0 |
|
532 |
|
533 set_line(lines,nline,header) |
|
534 set_line(lines,nline,table_header) |
|
535 ;----------------------------------------------- |
|
536 ;row of table |
|
537 |
|
538 do n = 0,nrow-1 |
|
539 set_line(lines,nline,row_header) |
|
540 |
|
541 txt0 = row_head(n) |
|
542 txt1 = text(n,0) |
|
543 txt2 = text(n,1) |
|
544 txt3 = text(n,2) |
|
545 txt4 = text(n,3) |
|
546 txt5 = text(n,4) |
|
547 txt6 = text(n,5) |
|
548 txt7 = text(n,6) |
|
549 txt8 = text(n,7) |
|
550 txt9 = text(n,8) |
|
551 |
|
552 set_line(lines,nline,"<th>"+txt0+"</th>") |
|
553 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
554 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
555 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
556 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
557 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
558 set_line(lines,nline,"<th>"+txt6+"</th>") |
|
559 set_line(lines,nline,"<th>"+txt7+"</th>") |
|
560 set_line(lines,nline,"<th>"+txt8+"</th>") |
|
561 set_line(lines,nline,"<th>"+txt9+"</th>") |
|
562 |
|
563 set_line(lines,nline,row_footer) |
|
564 end do |
|
565 ;----------------------------------------------- |
|
566 set_line(lines,nline,table_footer) |
|
567 set_line(lines,nline,footer) |
|
568 |
|
569 ; Now write to an HTML file. |
|
570 |
|
571 output_html = "table_per_biome.html" |
|
572 |
|
573 idx = ind(.not.ismissing(lines)) |
|
574 if(.not.any(ismissing(idx))) then |
|
575 asciiwrite(output_html,lines(idx)) |
|
576 else |
|
577 print ("error?") |
|
578 end if |
|
579 |
|
580 delete (idx) |
|
581 |
|
582 ;--------------------------------------------------- |
|
583 ; read model data, time series: |
|
584 |
|
585 fm = addfile (dirm+film7,"r") |
|
586 |
|
587 NPP3 = fm->NPP |
|
588 NEE3 = fm->NEE |
|
589 ;Fire = fm->COL_FIRE_CLOSS |
|
590 |
|
591 delete (fm) |
|
592 |
|
593 ; Units for these variables are: |
|
594 |
|
595 ;NPP3: g C/m^2/s |
|
596 ;NEE3: g C/m^2/s |
|
597 ;Fire: g C/m^2/s |
|
598 |
|
599 nsec_per_month = 60*60*24*30 |
|
600 |
|
601 ; change unit to g C/m^2/month |
|
602 |
|
603 NPP3 = NPP3 * nsec_per_month * conform(NPP3,landfrac,(/2,3/)) |
|
604 NEE3 = NEE3 * nsec_per_month * conform(NEE3,landfrac,(/2,3/)) |
|
605 ;Fire = Fire * nsec_per_month * conform(Fire,landfrac,(/2,3/)) |
|
606 |
|
607 ;data_n = 3 |
|
608 data_n = 2 |
|
609 |
|
610 dsizes = dimsizes(NPP3) |
|
611 nyear = dsizes(0) |
|
612 nmonth = dsizes(1) |
|
613 ntime = nyear * nmonth |
|
614 |
|
615 year_start = 1979 |
|
616 year_end = 2004 |
|
617 |
|
618 ;******************************************************************* |
|
619 ; Calculate "nice" bins for binning the data in equally spaced ranges |
|
620 ;******************************************************************** |
|
621 |
|
622 ; using model biome class |
|
623 nclass = nclass_mod |
|
624 |
|
625 range = fspan(0,nclass,nclass+1) |
|
626 |
|
627 ; print (range) |
|
628 ; Use this range information to grab all the values in a |
|
629 ; particular range, and then take an average. |
|
630 |
|
631 nx = dimsizes(range) - 1 |
|
632 |
|
633 ;============================== |
|
634 ; put data into bins |
|
635 ;============================== |
|
636 |
|
637 ; using observed biome class |
|
638 ; base = ndtooned(classob) |
|
639 ; using model biome class |
|
640 base = ndtooned(classmod) |
|
641 |
|
642 ; output |
|
643 |
|
644 yvalues = new((/ntime,data_n,nx/),float) |
|
645 |
|
646 ; Loop through each range, using base. |
|
647 |
|
648 do i=0,nx-1 |
|
649 |
|
650 if (i.ne.(nx-1)) then |
|
651 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) |
|
652 else |
|
653 idx = ind(base.ge.range(i)) |
|
654 end if |
|
655 |
|
656 do n = 0,data_n-1 |
|
657 |
|
658 t = -1 |
|
659 do m = 0,nyear-1 |
|
660 do k = 0,nmonth-1 |
|
661 |
|
662 t = t + 1 |
|
663 |
|
664 if (n.eq.0) then |
|
665 data = ndtooned(NPP3(m,k,:,:)) |
|
666 end if |
|
667 |
|
668 if (n.eq.1) then |
|
669 data = ndtooned(NEE3(m,k,:,:)) |
|
670 end if |
|
671 |
|
672 ; if (n.eq.2) then |
|
673 ; data = ndtooned(Fire(m,k,:,:)) |
|
674 ; end if |
|
675 |
|
676 ; Calculate average |
|
677 |
|
678 if (.not.any(ismissing(idx))) then |
|
679 yvalues(t,n,i) = sum(data(idx)*area_1d(idx)) |
|
680 else |
|
681 yvalues(t,n,i) = yvalues@_FillValue |
|
682 end if |
|
683 |
|
684 ;############################################################# |
|
685 ; using model biome class: |
|
686 ; set the following 4 classes to _FillValue: |
|
687 ; (3)Needleleaf Deciduous Boreal Tree, |
|
688 ; (8)Broadleaf Deciduous Boreal Tree, |
|
689 ; (9)Broadleaf Evergreen Shrub, |
|
690 ; (16)Wheat |
|
691 |
|
692 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then |
|
693 yvalues(t,n,i) = yvalues@_FillValue |
|
694 end if |
|
695 ;############################################################# |
|
696 |
|
697 end do |
|
698 end do |
|
699 |
|
700 delete(data) |
|
701 end do |
|
702 |
|
703 delete(idx) |
|
704 end do |
|
705 |
|
706 delete (base) |
|
707 delete (NPP3) |
|
708 delete (NEE3) |
|
709 ; delete (Fire) |
|
710 |
|
711 ;---------------------------------------------------------------- |
|
712 ; data for tseries plot |
|
713 |
|
714 yvalues_g = new((/ntime,data_n,n_biome/),float) |
|
715 |
|
716 yvalues_g@units = "TgC/month" |
|
717 |
|
718 ; change unit to Tg C/month |
|
719 ; change unit from g to Tg (Tera gram) |
|
720 factor_unit = 1.e-12 |
|
721 |
|
722 yvalues_g = yvalues(:,:,good) * factor_unit |
|
723 |
|
724 ;******************************************************************* |
|
725 ; general settings for line plot |
|
726 ;******************************************************************* |
|
727 |
|
728 ; res |
|
729 res = True |
|
730 res@xyDashPatterns = (/0/) ; make lines solid |
|
731 res@xyLineThicknesses = (/2.0/) ; make lines thicker |
|
732 res@xyLineColors = (/"blue"/) ; line color |
|
733 |
|
734 res@trXMinF = year_start |
|
735 res@trXMaxF = year_end + 1 |
|
736 |
|
737 res@vpHeightF = 0.4 ; change aspect ratio of plot |
|
738 ; res@vpWidthF = 0.8 |
|
739 res@vpWidthF = 0.75 |
|
740 |
|
741 ; res@gsnMaximize = True |
|
742 |
|
743 ;******************************************************************* |
|
744 ; (A) 1 component in each biome: monthly |
|
745 ;******************************************************************* |
|
746 |
|
747 ; component = (/"NPP","NEE","Fire"/) |
|
748 component = (/"NPP","NEE"/) |
|
749 |
|
750 ; for x-axis in xyplot |
|
751 |
|
752 timeI = new((/ntime/),integer) |
|
753 timeF = new((/ntime/),float) |
|
754 timeI = ispan(1,ntime,1) |
|
755 timeF = year_start + (timeI-1)/12. |
|
756 timeF@long_name = "year" |
|
757 |
|
758 plot_data = new((/ntime/),float) |
|
759 plot_data@long_name = "TgC/month" |
|
760 |
|
761 do n = 0, data_n-1 |
|
762 do m = 0, n_biome-1 |
|
763 |
|
764 plot_name = component(n)+"_monthly_biome_"+ m |
|
765 |
|
766 wks = gsn_open_wks (plot_type,plot_name) |
|
767 |
|
768 title = component(n)+ ": "+ row_head(m) |
|
769 res@tiMainString = title |
|
770 res@tiMainFontHeightF = 0.025 |
|
771 |
|
772 plot_data(:) = yvalues_g(:,n,m) |
|
773 |
|
774 plot=gsn_csm_xy(wks,timeF,plot_data,res) |
|
775 |
|
776 delete (wks) |
|
777 delete (plot) |
|
778 |
|
779 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
780 "rm "+plot_name+"."+plot_type) |
|
781 end do |
|
782 end do |
|
783 |
|
784 do n = 0, data_n-1 |
|
785 |
|
786 plot_name = component(n)+"_monthly_global" |
|
787 |
|
788 wks = gsn_open_wks (plot_type,plot_name) |
|
789 |
|
790 title = component(n)+ ": Global" |
|
791 res@tiMainString = title |
|
792 res@tiMainFontHeightF = 0.025 |
|
793 |
|
794 do k = 0,ntime-1 |
|
795 plot_data(k) = sum(yvalues_g(k,n,:)) |
|
796 end do |
|
797 |
|
798 plot=gsn_csm_xy(wks,timeF,plot_data,res) |
|
799 |
|
800 delete (wks) |
|
801 delete (plot) |
|
802 |
|
803 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
804 "rm "+plot_name+"."+plot_type) |
|
805 end do |
|
806 |
|
807 delete (plot_data) |
|
808 delete (timeI) |
|
809 delete (timeF) |
|
810 |
|
811 ;******************************************************************* |
|
812 ; (B) 1 component in each biome: annually |
|
813 ;******************************************************************* |
|
814 |
|
815 yvalues_a = new((/nyear,data_n,n_biome/),float) |
|
816 yvalues_g!0 = "time" |
|
817 yvalues_g!1 = "case" |
|
818 yvalues_g!2 = "record" |
|
819 |
|
820 yvalues_a = month_to_annual(yvalues_g,0) |
|
821 |
|
822 delete (yvalues_g) |
|
823 |
|
824 ; for x-axis in xyplot |
|
825 |
|
826 timeI = new((/nyear/),integer) |
|
827 timeF = new((/nyear/),float) |
|
828 timeI = ispan(1,nyear,1) |
|
829 timeF = year_start + (timeI-1) |
|
830 timeF@long_name = "year" |
|
831 |
|
832 plot_data = new((/nyear/),float) |
|
833 plot_data@long_name = "TgC/year" |
|
834 |
|
835 do n = 0, data_n-1 |
|
836 do m = 0, n_biome-1 |
|
837 |
|
838 plot_name = component(n)+"_annual_biome_"+ m |
|
839 |
|
840 wks = gsn_open_wks (plot_type,plot_name) |
|
841 |
|
842 title = component(n)+ ": "+ row_head(m) |
|
843 res@tiMainString = title |
|
844 res@tiMainFontHeightF = 0.025 |
|
845 |
|
846 plot_data(:) = yvalues_a(:,n,m) |
|
847 |
|
848 plot=gsn_csm_xy(wks,timeF,plot_data,res) |
|
849 |
|
850 delete (wks) |
|
851 delete (plot) |
|
852 |
|
853 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
854 "rm "+plot_name+"."+plot_type) |
|
855 end do |
|
856 end do |
|
857 |
|
858 do n = 0, data_n-1 |
|
859 |
|
860 plot_name = component(n)+"_annual_global" |
|
861 |
|
862 wks = gsn_open_wks (plot_type,plot_name) |
|
863 |
|
864 title = component(n)+ ": Global" |
|
865 res@tiMainString = title |
|
866 res@tiMainFontHeightF = 0.025 |
|
867 |
|
868 do k = 0,nyear-1 |
|
869 plot_data(k) = sum(yvalues_a(k,n,:)) |
|
870 end do |
|
871 |
|
872 plot=gsn_csm_xy(wks,timeF,plot_data,res) |
|
873 |
|
874 delete (wks) |
|
875 delete (plot) |
|
876 |
|
877 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
878 "rm "+plot_name+"."+plot_type) |
|
879 end do |
|
880 |
|
881 ;**************************************** |
|
882 ; output plot and html |
|
883 ;**************************************** |
|
884 output_dir = model_name+"/carbon_sink" |
|
885 |
|
886 system("mv *.png *.html " + output_dir) |
|
887 ;**************************************** |
|
888 |
|
889 end |
|
890 |