all/04.biomass.ncl
changeset 1 4be95183fbcd
parent 0 0c6405ab2ff4
equal deleted inserted replaced
0:7e30d378d3f4 1:d253bb5fd427
    42   area0     = area0 * landfrac
    42   area0     = area0 * landfrac
    43 
    43 
    44 ;--------------------------------------------
    44 ;--------------------------------------------
    45 ; read model data
    45 ; read model data
    46 
    46 
    47  fm   = addfile (dirm+film1,"r")
    47  fm   = addfile (dirm+film11,"r")
    48 
    48 
    49  if (BGC .eq. "cn") then 
    49  if (BGC .eq. "cn") then 
    50     data1  = fm->LIVESTEMC
    50     data1  = fm->LIVESTEMC
    51     data2  = fm->DEADSTEMC
    51     data2  = fm->DEADSTEMC
    52     data3  = fm->LEAFC
    52     data3  = fm->LEAFC
    81   delete (fm)
    81   delete (fm)
    82 
    82 
    83 ;------------------------------------------------
    83 ;------------------------------------------------
    84 ; read ob data
    84 ; read ob data
    85 
    85 
    86  ob_name = "LC15_Amazon_Biomass"
    86  ob_name = "LC15 Amazon Biomass"
    87 
    87 
    88  dir_b = diro + "biomass/"
    88  dir_b = diro + "biomass/"
    89  fil_b = "amazon_biomass_"+model_grid+".nc"
    89  fil_b = "amazon_biomass_"+model_grid+".nc"
    90  fo   = addfile (dir_b+fil_b,"r")
    90  fo   = addfile (dir_b+fil_b,"r")
    91  
    91  
    97 
    97 
    98 ;************************************************
    98 ;************************************************
    99 ; Units for these variables are:
    99 ; Units for these variables are:
   100 ; dataob   : MgC/ha
   100 ; dataob   : MgC/ha
   101 ; datamod0 : gC/m2
   101 ; datamod0 : gC/m2
   102 ; We want to convert these to KgC/m2
   102 ; We want to convert these to kgC/m2
   103 ; ha = 100m*100m = 10,000 m2
   103 ; ha = 100m*100m = 10,000 m2
   104 ; MgC/ha*1000/10,000 = KgC/m2
   104 ; MgC/ha*1000/10,000 = kgC/m2
   105 
   105 
   106   factor_aboveground = 0.5
   106   factor_aboveground = 0.5
   107   factor_unit_ob     = 0.1
   107   factor_unit_ob     = 0.1
   108   factor_unit_mod    = 0.001         
   108   factor_unit_mod    = 0.001         
   109 
   109 
   110   dataob   = dataob * factor_aboveground * factor_unit_ob
   110   dataob   = dataob * factor_aboveground * factor_unit_ob
   111   datamod0 = datamod0 * factor_unit_mod 
   111   datamod0 = datamod0 * factor_unit_mod 
   112 
   112 
   113   dataob@units      = "KgC/m2"
   113   dataob@units      = "kg C m~S~-2~N~"
   114   datamod0@units    = "KgC/m2"
   114   datamod0@units    = "kg C m~S~-2~N~"
   115 
   115 
   116   dataob@long_name      = "Amazon Biomass"
   116   dataob@long_name      = "Amazon Biomass"
   117   datamod0@long_name    = "Amazon Biomass"
   117   datamod0@long_name    = "Amazon Biomass"
   118 ;********************************************************
   118 ;********************************************************
   119 ; get subset of datamod0 that match dataob
   119 ; get subset of datamod0 that match dataob
   137 
   137 
   138   mask_amazon@units = ""
   138   mask_amazon@units = ""
   139 ;********************************************************
   139 ;********************************************************
   140 ; sum over amazom_mask area:
   140 ; sum over amazom_mask area:
   141 
   141 
   142 ; Peta g = 1.e15 g = 1.e12 Kg
   142 ; Peta g = 1.e15 g = 1.e12 kg
   143   factor_unit = 1.e-12
   143   factor_unit = 1.e-12
   144 
   144 
   145 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
   145 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
   146 
   146 
   147   Sum_area = sum(area*mask_amazon)*factor_unit
   147   Sum_area = sum(area*mask_amazon)*factor_unit
   202 
   202 
   203   delete (wks)
   203   delete (wks)
   204 
   204 
   205   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   205   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   206          "rm "+plot_name+"."+plot_type)
   206          "rm "+plot_name+"."+plot_type)
       
   207   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   207 
   208 
   208 ;------------------------------------------------------------------------
   209 ;------------------------------------------------------------------------
   209 ; contour ob
   210 ; contour ob
   210 
   211 
   211   resg@cnMinLevelValF       = 0.              ; Min level
   212   resg@cnMinLevelValF       = 0.              ; Min level
   223 
   224 
   224   delete (wks)
   225   delete (wks)
   225 
   226 
   226   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   227   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   227          "rm "+plot_name+"."+plot_type)
   228          "rm "+plot_name+"."+plot_type)
       
   229   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   228 
   230 
   229 ;------------------------------------------------------------------------
   231 ;------------------------------------------------------------------------
   230 ; contour model
   232 ; contour model
   231 
   233 
   232   resg@cnMinLevelValF       = 0.              ; Min level
   234   resg@cnMinLevelValF       = 0.              ; Min level
   233   resg@cnMaxLevelValF       = 30.             ; Max level
   235   resg@cnMaxLevelValF       = 30.             ; Max level
   234   resg@cnLevelSpacingF      = 2.              ; interval
   236   resg@cnLevelSpacingF      = 2.              ; interval
   235 
   237 
   236   plot_name = "global_model"
   238   plot_name = "global_model"
   237   title     = "Model "+ model_name 
   239   title     = "Model "+ model_name + " (2000)"
   238   resg@tiMainString  = title
   240   resg@tiMainString  = title
   239 
   241 
   240   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
   242   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
   241   gsn_define_colormap(wks,"gui_default")     ; choose colormap
   243   gsn_define_colormap(wks,"gui_default")     ; choose colormap
   242 
   244 
   244 
   246 
   245   delete (wks)
   247   delete (wks)
   246 
   248 
   247   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   249   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   248          "rm "+plot_name+"."+plot_type)
   250          "rm "+plot_name+"."+plot_type)
       
   251   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   249 
   252 
   250 ;------------------------------------------------------------------------
   253 ;------------------------------------------------------------------------
   251 ; contour model vs ob
   254 ; contour model vs ob
   252 
   255 
   253   plot_name = "global_model_vs_ob"
   256   plot_name = "global_model_vs_ob"
   271   ug = uu(good)
   274   ug = uu(good)
   272   vg = vv(good)
   275   vg = vv(good)
   273 
   276 
   274   ccrG = esccr(ug,vg,0)
   277   ccrG = esccr(ug,vg,0)
   275 
   278 
   276   score_max = 5.0
   279   score_max = 10.0
   277 
   280 
   278 ; Miomass = (ccrG*ccrG)* score_max
   281 ; Miomass = (ccrG*ccrG)* score_max
   279 ; new eq
   282 ; new eq
   280   bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
   283   bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
       
   284   print ("bias=" + bias + "dimsizes(ug)=" + dimsizes(ug))
       
   285   print ("other bias=" + sum(abs(ug-vg)/(ug+vg)))
   281   Mbiomass  = (1. - (bias/dimsizes(ug)))*score_max
   286   Mbiomass  = (1. - (bias/dimsizes(ug)))*score_max
       
   287   print ("Mbiomass=" + Mbiomass)
   282   M_biomass = sprintf("%.2f", Mbiomass)
   288   M_biomass = sprintf("%.2f", Mbiomass)
       
   289   print ("M_biomass=" + M_biomass)
   283 
   290 
   284   if (isvar("compare")) then
   291   if (isvar("compare")) then
   285      system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \ 
   292      system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \ 
   286             "mv -f "+html_new2+" "+html_name2)
   293             "mv -f "+html_new2+" "+html_name2)
   287   end if
   294   end if
   307 
   314 
   308   plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)       
   315   plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)       
   309 
   316 
   310 ;(b) model
   317 ;(b) model
   311 
   318 
   312   title     = "Model "+ model_name
   319   title     = "Model "+ model_name + " (2000)"
   313   resg@tiMainString  = title
   320   resg@tiMainString  = title
   314 
   321 
   315   plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) 
   322   plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) 
   316 
   323 
   317 ;(c) model-ob
   324 ;(c) model-ob
   318 
   325 
   319   zz = datamod
   326   zz = datamod
   320   zz = datamod - dataob
   327   zz = datamod - dataob
   321   title = "Model_"+model_name+" - Observed"
   328   title = "Model "+model_name+" (2000) - Observed"
   322 
   329 
   323   resg@cnMinLevelValF  = -10.          ; Min level
   330   resg@cnMinLevelValF  = -10.          ; Min level
   324   resg@cnMaxLevelValF  =  10.          ; Max level
   331   resg@cnMaxLevelValF  =  10.          ; Max level
   325   resg@cnLevelSpacingF =  2.           ; interval
   332   resg@cnLevelSpacingF =  2.           ; interval
   326   resg@tiMainString    = title
   333   resg@tiMainString    = title
   338   delete (plot)
   345   delete (plot)
   339   delete (zz)
   346   delete (zz)
   340 
   347 
   341   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   348   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   342          "rm "+plot_name+"."+plot_type)
   349          "rm "+plot_name+"."+plot_type)
       
   350   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   343 
   351 
   344   resg@gsnFrame             = True            ; draw plot 
   352   resg@gsnFrame             = True            ; draw plot 
   345   resg@gsnDraw              = True            ; advance frame
   353   resg@gsnDraw              = True            ; advance frame
   346 ;------------------------------------------------------------------------
   354 ;------------------------------------------------------------------------
   347 ; contour ob : masked
   355 ; contour ob : masked
   361 
   369 
   362   gRes  = True
   370   gRes  = True
   363   gRes@txFontHeightF = 0.02
   371   gRes@txFontHeightF = 0.02
   364   gRes@txAngleF = 0
   372   gRes@txAngleF = 0
   365 
   373 
   366   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" Kg C/m2)"
   374   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" kg C m~S~-2~N~)"
   367 
   375 
   368   gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
   376   gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
   369 ;-----------------------------------------
   377 ;-----------------------------------------
   370   zo = dataob
   378   zo = dataob
   371   zo = dataob*mask_amazon
   379   zo = dataob*mask_amazon
   375   delete (wks)
   383   delete (wks)
   376   delete (plot)
   384   delete (plot)
   377 
   385 
   378   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   386   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   379          "rm "+plot_name+"."+plot_type)
   387          "rm "+plot_name+"."+plot_type)
       
   388   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   380 
   389 
   381 ;------------------------------------------------------------------------
   390 ;------------------------------------------------------------------------
   382 ; contour model: masked
   391 ; contour model: masked
   383 
   392 
   384   resg@cnMinLevelValF       = 0.              ; Min level
   393   resg@cnMinLevelValF       = 0.              ; Min level
   385   resg@cnMaxLevelValF       = 30.             ; Max level
   394   resg@cnMaxLevelValF       = 30.             ; Max level
   386   resg@cnLevelSpacingF      = 2.              ; interval
   395   resg@cnLevelSpacingF      = 2.              ; interval
   387 
   396 
   388   plot_name = "global_mask_model"
   397   plot_name = "global_mask_model"
   389   title     = "Model "+ model_name 
   398   title     = "Model "+ model_name + " (2000)"
   390   resg@tiMainString  = title
   399   resg@tiMainString  = title
   391 
   400 
   392   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
   401   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
   393   gsn_define_colormap(wks,"gui_default")     ; choose colormap
   402   gsn_define_colormap(wks,"gui_default")     ; choose colormap
   394 ;-----------------------------------------
   403 ;-----------------------------------------
   396 
   405 
   397   gRes  = True
   406   gRes  = True
   398   gRes@txFontHeightF = 0.02
   407   gRes@txFontHeightF = 0.02
   399   gRes@txAngleF = 0
   408   gRes@txAngleF = 0
   400 
   409 
   401   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" Kg C/m2)"
   410   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" kg C m~S~-2~N~)"
   402 
   411 
   403   gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
   412   gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
   404 ;-----------------------------------------
   413 ;-----------------------------------------
   405   zm = datamod
   414   zm = datamod
   406   zm = datamod*mask_amazon
   415   zm = datamod*mask_amazon
   410   delete (wks)
   419   delete (wks)
   411   delete (plot)
   420   delete (plot)
   412 
   421 
   413   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   422   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   414          "rm "+plot_name+"."+plot_type)
   423          "rm "+plot_name+"."+plot_type)
       
   424   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   415 
   425 
   416 ;------------------------------------------------------------------------
   426 ;------------------------------------------------------------------------
   417 ; contour model vs ob: masked 
   427 ; contour model vs ob: masked 
   418 
   428 
   419   plot_name = "global_mask_vs_ob"
   429   plot_name = "global_mask_vs_ob"
   432   delete (uu)
   442   delete (uu)
   433   delete (vv)
   443   delete (vv)
   434   delete (ug)
   444   delete (ug)
   435   delete (vg)  
   445   delete (vg)  
   436 
   446 
   437   score_max = 5.
   447   score_max = 0.
   438 
   448 
   439   uu = ndtooned(zm)
   449   uu = ndtooned(zm)
   440   vv = ndtooned(zo)
   450   vv = ndtooned(zo)
   441 
   451 
   442   good = ind((uu .gt. 0.) .and. (vv .gt. 0.))
   452   good = ind((uu .gt. 0.) .and. (vv .gt. 0.))
   486 
   496 
   487   plot(0) = gsn_csm_contour_map_ce(wks,zo,resg)       
   497   plot(0) = gsn_csm_contour_map_ce(wks,zo,resg)       
   488 
   498 
   489 ;(b) model
   499 ;(b) model
   490 
   500 
   491   title     = "Model "+ model_name
   501   title     = "Model "+ model_name + " (2000)"
   492   resg@tiMainString  = title
   502   resg@tiMainString  = title
   493 
   503 
   494   plot(1) = gsn_csm_contour_map_ce(wks,zm,resg) 
   504   plot(1) = gsn_csm_contour_map_ce(wks,zm,resg) 
   495 
   505 
   496 ;(c) model-ob
   506 ;(c) model-ob
   497 
   507 
   498   zz = zo
   508   zz = zo
   499   zz = zm - zo
   509   zz = zm - zo
   500   title = "Model_"+model_name+" - Observed"
   510   title = "Model "+model_name+" (2000) - Observed"
   501 
   511 
   502   resg@cnMinLevelValF  = -10.          ; Min level
   512   resg@cnMinLevelValF  = -10.          ; Min level
   503   resg@cnMaxLevelValF  =  10.          ; Max level
   513   resg@cnMaxLevelValF  =  10.          ; Max level
   504   resg@cnLevelSpacingF =  2.           ; interval
   514   resg@cnLevelSpacingF =  2.           ; interval
   505   resg@tiMainString    = title
   515   resg@tiMainString    = title
   516   delete (wks)
   526   delete (wks)
   517   delete (plot)
   527   delete (plot)
   518 
   528 
   519   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   529   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
   520          "rm "+plot_name+"."+plot_type)
   530          "rm "+plot_name+"."+plot_type)
       
   531   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   521 
   532 
   522 ;***************************************************************************
   533 ;***************************************************************************
   523 ; add total score and write to file
   534 ; add total score and write to file
   524 ;***************************************************************************
   535 ;***************************************************************************
   525   M_total = Mbiomass + Mbiomass2
   536   M_total = Mbiomass + Mbiomass2