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