|
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 begin |
|
7 |
|
8 plot_type = "ps" |
|
9 plot_type_new = "png" |
|
10 |
|
11 ;------------------------------------------------------ |
|
12 ; edit table.html of current model for movel1_vs_model2 |
|
13 |
|
14 if (isvar("compare")) then |
|
15 html_name2 = compare+"/table.html" |
|
16 html_new2 = html_name2 +".new" |
|
17 end if |
|
18 |
|
19 ;------------------------------------------------------ |
|
20 ; edit table.html for current model |
|
21 |
|
22 html_name = model_name+"/table.html" |
|
23 html_new = html_name +".new" |
|
24 |
|
25 ;--------------------------------- |
|
26 ; get model data: landfrac and area |
|
27 |
|
28 film_l = "lnd_"+ model_grid + ".nc" |
|
29 fm_l = addfile (dirs+film_l,"r") |
|
30 |
|
31 landfrac = fm_l->landfrac |
|
32 area0 = fm_l->area |
|
33 |
|
34 delete (fm_l) |
|
35 |
|
36 ; change area from km**2 to m**2 |
|
37 area0 = area0 * 1.e6 |
|
38 |
|
39 ;----------------------------- |
|
40 ; take into account landfrac |
|
41 |
|
42 area0 = area0 * landfrac |
|
43 |
|
44 ;-------------------------------------------- |
|
45 ; read model data |
|
46 |
|
47 fm = addfile (dirm+film1,"r") |
|
48 |
|
49 if (BGC .eq. "cn") then |
|
50 data1 = fm->LIVESTEMC |
|
51 data2 = fm->DEADSTEMC |
|
52 data3 = fm->LEAFC |
|
53 datamod0 = data1(0,:,:) |
|
54 datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:) |
|
55 end if |
|
56 |
|
57 if (BGC .eq. "casa") then |
|
58 factor_WOODC = 0.7 |
|
59 data1 = fm->WOODC |
|
60 data2 = fm->LEAFC |
|
61 datamod0 = data1(0,:,:) |
|
62 datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:) |
|
63 end if |
|
64 |
|
65 ; unit: gC/m2 |
|
66 |
|
67 xm = fm->lon |
|
68 ym = fm->lat |
|
69 |
|
70 delete (fm) |
|
71 |
|
72 ;------------------------------------------------ |
|
73 ; read amazon mask data |
|
74 |
|
75 dir_m = diro + "biomass/" |
|
76 fil_m = "amazon_mask_"+ model_grid + ".nc" |
|
77 fm = addfile (dir_m+fil_m,"r") |
|
78 |
|
79 mask_amazon0 = fm->mask_amazon |
|
80 |
|
81 delete (fm) |
|
82 |
|
83 ;------------------------------------------------ |
|
84 ; read ob data |
|
85 |
|
86 ob_name = "LC15_Amazon_Biomass" |
|
87 |
|
88 dir_b = diro + "biomass/" |
|
89 fil_b = "amazon_biomass_"+model_grid+".nc" |
|
90 fo = addfile (dir_b+fil_b,"r") |
|
91 |
|
92 dataob = fo->BIOMASS |
|
93 xo = fo->lon |
|
94 yo = fo->lat |
|
95 |
|
96 delete (fo) |
|
97 |
|
98 ;************************************************ |
|
99 ; Units for these variables are: |
|
100 ; dataob : MgC/ha |
|
101 ; datamod0 : gC/m2 |
|
102 ; We want to convert these to KgC/m2 |
|
103 ; ha = 100m*100m = 10,000 m2 |
|
104 ; MgC/ha*1000/10,000 = KgC/m2 |
|
105 |
|
106 factor_aboveground = 0.5 |
|
107 factor_unit_ob = 0.1 |
|
108 factor_unit_mod = 0.001 |
|
109 |
|
110 dataob = dataob * factor_aboveground * factor_unit_ob |
|
111 datamod0 = datamod0 * factor_unit_mod |
|
112 |
|
113 dataob@units = "KgC/m2" |
|
114 datamod0@units = "KgC/m2" |
|
115 |
|
116 dataob@long_name = "Amazon Biomass" |
|
117 datamod0@long_name = "Amazon Biomass" |
|
118 ;******************************************************** |
|
119 ; get subset of datamod0 that match dataob |
|
120 |
|
121 nlon = dimsizes(xo) |
|
122 nlat = dimsizes(yo) |
|
123 |
|
124 ind_lonL = ind(xm .eq. xo(0)) |
|
125 ind_lonR = ind(xm .eq. xo(nlon-1)) |
|
126 ind_latS = ind(ym .eq. yo(0)) |
|
127 ind_latN = ind(ym .eq. yo(nlat-1)) |
|
128 |
|
129 datamod = dataob |
|
130 datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR) |
|
131 |
|
132 area = dataob |
|
133 area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR) |
|
134 |
|
135 mask_amazon = dataob |
|
136 mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR) |
|
137 |
|
138 mask_amazon@units = "" |
|
139 ;******************************************************** |
|
140 ; sum over amazom_mask area: |
|
141 |
|
142 ; Peta g = 1.e15 g = 1.e12 Kg |
|
143 factor_unit = 1.e-12 |
|
144 |
|
145 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.) |
|
146 |
|
147 Sum_area = sum(area*mask_amazon)*factor_unit |
|
148 |
|
149 Sum_ob = sum(dataob*area*mask_amazon)*factor_unit |
|
150 Sum_mod = sum(datamod*area*mask_amazon)*factor_unit |
|
151 |
|
152 avg_ob = Sum_ob /Sum_area |
|
153 avg_mod = Sum_mod/Sum_area |
|
154 |
|
155 Sum_biomass_ob = sprintf("%.2f",Sum_ob ) |
|
156 Sum_biomass_mod = sprintf("%.2f",Sum_mod) |
|
157 |
|
158 ;---------------------------------------------------------------------- |
|
159 ; contour plot res |
|
160 |
|
161 resg = True ; Use plot options |
|
162 resg@cnFillOn = True ; Turn on color fill |
|
163 resg@gsnSpreadColors = True ; use full colormap |
|
164 resg@cnLinesOn = False ; Turn off contourn lines |
|
165 resg@mpFillOn = False ; Turn off map fill |
|
166 resg@gsnAddCyclic = False |
|
167 |
|
168 resg@gsnSpreadColors = True ; use full colormap |
|
169 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
170 |
|
171 resg@mpMinLatF = -21.1 ; range to zoom in on |
|
172 resg@mpMaxLatF = 13.8 |
|
173 resg@mpMinLonF = 277.28 |
|
174 resg@mpMaxLonF = 326.43 |
|
175 ;------------------------------------------------------------------------ |
|
176 ; mask plot |
|
177 |
|
178 plot_name = "mask_ob" |
|
179 |
|
180 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
181 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
182 |
|
183 ;----------------------------------------- |
|
184 ; plot area sum |
|
185 |
|
186 gRes = True |
|
187 gRes@txFontHeightF = 0.02 |
|
188 ; gRes@txAngleF = 90 |
|
189 |
|
190 area_sum_text = "(mask area = "+sprintf("%.2f", Sum_area)+"e12 m2)" |
|
191 |
|
192 gsn_text_ndc(wks,area_sum_text,0.50,0.80,gRes) |
|
193 ;----------------------------------------- |
|
194 |
|
195 resg@cnMinLevelValF = 0. ; Min level |
|
196 resg@cnMaxLevelValF = 1. ; Max level |
|
197 resg@cnLevelSpacingF = 0.1 ; interval |
|
198 |
|
199 resg@tiMainString = "Amazon Mask: grid = "+ model_grid |
|
200 |
|
201 plot = gsn_csm_contour_map_ce(wks,mask_amazon,resg) |
|
202 |
|
203 delete (wks) |
|
204 |
|
205 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
206 "rm "+plot_name+"."+plot_type) |
|
207 |
|
208 ;------------------------------------------------------------------------ |
|
209 ; contour ob |
|
210 |
|
211 resg@cnMinLevelValF = 0. ; Min level |
|
212 resg@cnMaxLevelValF = 30. ; Max level |
|
213 resg@cnLevelSpacingF = 2. ; interval |
|
214 |
|
215 plot_name = "global_ob" |
|
216 title = ob_name+" "+ model_grid |
|
217 resg@tiMainString = title |
|
218 |
|
219 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
220 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
221 |
|
222 plot = gsn_csm_contour_map_ce(wks,dataob,resg) |
|
223 |
|
224 delete (wks) |
|
225 |
|
226 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
227 "rm "+plot_name+"."+plot_type) |
|
228 |
|
229 ;------------------------------------------------------------------------ |
|
230 ; contour model |
|
231 |
|
232 resg@cnMinLevelValF = 0. ; Min level |
|
233 resg@cnMaxLevelValF = 30. ; Max level |
|
234 resg@cnLevelSpacingF = 2. ; interval |
|
235 |
|
236 plot_name = "global_model" |
|
237 title = "Model "+ model_name |
|
238 resg@tiMainString = title |
|
239 |
|
240 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
241 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
242 |
|
243 plot = gsn_csm_contour_map_ce(wks,datamod,resg) |
|
244 |
|
245 delete (wks) |
|
246 |
|
247 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
248 "rm "+plot_name+"."+plot_type) |
|
249 |
|
250 ;------------------------------------------------------------------------ |
|
251 ; contour model vs ob |
|
252 |
|
253 plot_name = "global_model_vs_ob" |
|
254 |
|
255 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
256 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
257 |
|
258 delete (plot) |
|
259 plot=new(3,graphic) ; create graphic array |
|
260 |
|
261 resg@gsnFrame = False ; Do not draw plot |
|
262 resg@gsnDraw = False ; Do not advance frame |
|
263 |
|
264 ;(d) compute correlation coef and M score |
|
265 |
|
266 uu = ndtooned(datamod) |
|
267 vv = ndtooned(dataob) |
|
268 |
|
269 good = ind(.not.ismissing(uu) .and. .not.ismissing(vv)) |
|
270 |
|
271 ug = uu(good) |
|
272 vg = vv(good) |
|
273 |
|
274 ccrG = esccr(ug,vg,0) |
|
275 |
|
276 score_max = 5.0 |
|
277 |
|
278 ; Miomass = (ccrG*ccrG)* score_max |
|
279 ; new eq |
|
280 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg))) |
|
281 Mbiomass = (1. - (bias/dimsizes(ug)))*score_max |
|
282 M_biomass = sprintf("%.2f", Mbiomass) |
|
283 |
|
284 if (isvar("compare")) then |
|
285 system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
286 "mv -f "+html_new2+" "+html_name2) |
|
287 end if |
|
288 |
|
289 system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \ |
|
290 "mv -f "+html_new+" "+html_name) |
|
291 |
|
292 ; plot correlation coef |
|
293 |
|
294 gRes = True |
|
295 gRes@txFontHeightF = 0.02 |
|
296 gRes@txAngleF = 90 |
|
297 |
|
298 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")" |
|
299 |
|
300 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
301 ;-------------------------------------------------------------------- |
|
302 |
|
303 ;(a) ob |
|
304 |
|
305 title = ob_name+" "+ model_grid |
|
306 resg@tiMainString = title |
|
307 |
|
308 plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg) |
|
309 |
|
310 ;(b) model |
|
311 |
|
312 title = "Model "+ model_name |
|
313 resg@tiMainString = title |
|
314 |
|
315 plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) |
|
316 |
|
317 ;(c) model-ob |
|
318 |
|
319 zz = datamod |
|
320 zz = datamod - dataob |
|
321 title = "Model_"+model_name+" - Observed" |
|
322 |
|
323 resg@cnMinLevelValF = -10. ; Min level |
|
324 resg@cnMaxLevelValF = 10. ; Max level |
|
325 resg@cnLevelSpacingF = 2. ; interval |
|
326 resg@tiMainString = title |
|
327 |
|
328 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
329 |
|
330 pres = True ; panel plot mods desired |
|
331 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
332 ; indiv. plots in panel |
|
333 pres@gsnMaximize = True ; fill the page |
|
334 |
|
335 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
336 |
|
337 delete (wks) |
|
338 delete (plot) |
|
339 delete (zz) |
|
340 |
|
341 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
342 "rm "+plot_name+"."+plot_type) |
|
343 |
|
344 resg@gsnFrame = True ; draw plot |
|
345 resg@gsnDraw = True ; advance frame |
|
346 ;------------------------------------------------------------------------ |
|
347 ; contour ob : masked |
|
348 |
|
349 resg@cnMinLevelValF = 0. ; Min level |
|
350 resg@cnMaxLevelValF = 30. ; Max level |
|
351 resg@cnLevelSpacingF = 2. ; interval |
|
352 |
|
353 plot_name = "global_mask_ob" |
|
354 title = ob_name+" "+ model_grid |
|
355 resg@tiMainString = title |
|
356 |
|
357 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
358 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
359 ;----------------------------------------- |
|
360 ; plot average over mask region |
|
361 |
|
362 gRes = True |
|
363 gRes@txFontHeightF = 0.02 |
|
364 gRes@txAngleF = 0 |
|
365 |
|
366 area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" Kg C/m2)" |
|
367 |
|
368 gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes) |
|
369 ;----------------------------------------- |
|
370 zo = dataob |
|
371 zo = dataob*mask_amazon |
|
372 zo = where((mask_amazon .le. 0.01), zo@_FillValue, zo) |
|
373 plot = gsn_csm_contour_map_ce(wks,zo,resg) |
|
374 |
|
375 delete (wks) |
|
376 delete (plot) |
|
377 |
|
378 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
379 "rm "+plot_name+"."+plot_type) |
|
380 |
|
381 ;------------------------------------------------------------------------ |
|
382 ; contour model: masked |
|
383 |
|
384 resg@cnMinLevelValF = 0. ; Min level |
|
385 resg@cnMaxLevelValF = 30. ; Max level |
|
386 resg@cnLevelSpacingF = 2. ; interval |
|
387 |
|
388 plot_name = "global_mask_model" |
|
389 title = "Model "+ model_name |
|
390 resg@tiMainString = title |
|
391 |
|
392 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
393 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
394 ;----------------------------------------- |
|
395 ; plot average over mask region |
|
396 |
|
397 gRes = True |
|
398 gRes@txFontHeightF = 0.02 |
|
399 gRes@txAngleF = 0 |
|
400 |
|
401 area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" Kg C/m2)" |
|
402 |
|
403 gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes) |
|
404 ;----------------------------------------- |
|
405 zm = datamod |
|
406 zm = datamod*mask_amazon |
|
407 zm = where((mask_amazon .le. 0.01), zm@_FillValue, zm) |
|
408 plot = gsn_csm_contour_map_ce(wks,zm,resg) |
|
409 |
|
410 delete (wks) |
|
411 delete (plot) |
|
412 |
|
413 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
414 "rm "+plot_name+"."+plot_type) |
|
415 |
|
416 ;------------------------------------------------------------------------ |
|
417 ; contour model vs ob: masked |
|
418 |
|
419 plot_name = "global_mask_vs_ob" |
|
420 |
|
421 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
422 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
423 |
|
424 plot=new(3,graphic) ; create graphic array |
|
425 |
|
426 resg@gsnFrame = False ; Do not draw plot |
|
427 resg@gsnDraw = False ; Do not advance frame |
|
428 |
|
429 ;(d) compute correlation coef and M score |
|
430 |
|
431 delete (good) |
|
432 delete (uu) |
|
433 delete (vv) |
|
434 delete (ug) |
|
435 delete (vg) |
|
436 |
|
437 score_max = 5. |
|
438 |
|
439 uu = ndtooned(zm) |
|
440 vv = ndtooned(zo) |
|
441 |
|
442 good = ind((uu .gt. 0.) .and. (vv .gt. 0.)) |
|
443 |
|
444 ug = uu(good) |
|
445 vg = vv(good) |
|
446 |
|
447 ccrG = esccr(ug,vg,0) |
|
448 |
|
449 ; Miomass = (ccrG*ccrG)*score_max |
|
450 ; new eq |
|
451 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg))) |
|
452 Mbiomass2 = (1. - (bias/dimsizes(ug)))*score_max |
|
453 M_biomask = sprintf("%.2f", Mbiomass2) |
|
454 |
|
455 if (isvar("compare")) then |
|
456 system("sed -e '1,/M_biomask/s/M_biomask/"+M_biomask+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
457 "mv -f "+html_new2+" "+html_name2) |
|
458 system("sed -e '1,/Sum_biomass_ob/s/Sum_biomass_ob/"+Sum_biomass_ob+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
459 "mv -f "+html_new2+" "+html_name2) |
|
460 system("sed -e '1,/Sum_biomass_mod/s/Sum_biomass_mod/"+Sum_biomass_mod+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
461 "mv -f "+html_new2+" "+html_name2) |
|
462 end if |
|
463 |
|
464 system("sed s#M_biomask#"+M_biomask+"# "+html_name+" > "+html_new+";"+ \ |
|
465 "mv -f "+html_new+" "+html_name) |
|
466 system("sed s#Sum_biomass_ob#"+Sum_biomass_ob+"# "+html_name+" > "+html_new+";"+ \ |
|
467 "mv -f "+html_new+" "+html_name) |
|
468 system("sed s#Sum_biomass_mod#"+Sum_biomass_mod+"# "+html_name+" > "+html_new+";"+ \ |
|
469 "mv -f "+html_new+" "+html_name) |
|
470 ;-------------------------------------------------------------------- |
|
471 ; plot correlation coef |
|
472 |
|
473 gRes = True |
|
474 gRes@txFontHeightF = 0.02 |
|
475 gRes@txAngleF = 90 |
|
476 |
|
477 correlation_text = "(correlation coef = "+sprintf("%.2f", ccrG)+")" |
|
478 |
|
479 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
480 ;-------------------------------------------------------------------- |
|
481 |
|
482 ;(a) ob |
|
483 |
|
484 title = ob_name+" "+ model_grid |
|
485 resg@tiMainString = title |
|
486 |
|
487 plot(0) = gsn_csm_contour_map_ce(wks,zo,resg) |
|
488 |
|
489 ;(b) model |
|
490 |
|
491 title = "Model "+ model_name |
|
492 resg@tiMainString = title |
|
493 |
|
494 plot(1) = gsn_csm_contour_map_ce(wks,zm,resg) |
|
495 |
|
496 ;(c) model-ob |
|
497 |
|
498 zz = zo |
|
499 zz = zm - zo |
|
500 title = "Model_"+model_name+" - Observed" |
|
501 |
|
502 resg@cnMinLevelValF = -10. ; Min level |
|
503 resg@cnMaxLevelValF = 10. ; Max level |
|
504 resg@cnLevelSpacingF = 2. ; interval |
|
505 resg@tiMainString = title |
|
506 |
|
507 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
508 |
|
509 pres = True ; panel plot mods desired |
|
510 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
511 ; indiv. plots in panel |
|
512 pres@gsnMaximize = True ; fill the page |
|
513 |
|
514 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
515 |
|
516 delete (wks) |
|
517 delete (plot) |
|
518 |
|
519 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
520 "rm "+plot_name+"."+plot_type) |
|
521 |
|
522 ;*************************************************************************** |
|
523 ; add total score and write to file |
|
524 ;*************************************************************************** |
|
525 M_total = Mbiomass + Mbiomass2 |
|
526 |
|
527 asciiwrite("M_save.biomass", M_total) |
|
528 |
|
529 ;*************************************************************************** |
|
530 ; output plots |
|
531 ;*************************************************************************** |
|
532 output_dir = model_name+"/biomass" |
|
533 |
|
534 system("mv *.png " + output_dir) |
|
535 ;*************************************************************************** |
|
536 end |