1.1 --- a/all/04.biomass.ncl Mon Jan 26 22:08:20 2009 -0500
1.2 +++ b/all/04.biomass.ncl Thu Mar 26 14:02:21 2009 -0400
1.3 @@ -44,7 +44,7 @@
1.4 ;--------------------------------------------
1.5 ; read model data
1.6
1.7 - fm = addfile (dirm+film1,"r")
1.8 + fm = addfile (dirm+film11,"r")
1.9
1.10 if (BGC .eq. "cn") then
1.11 data1 = fm->LIVESTEMC
1.12 @@ -83,7 +83,7 @@
1.13 ;------------------------------------------------
1.14 ; read ob data
1.15
1.16 - ob_name = "LC15_Amazon_Biomass"
1.17 + ob_name = "LC15 Amazon Biomass"
1.18
1.19 dir_b = diro + "biomass/"
1.20 fil_b = "amazon_biomass_"+model_grid+".nc"
1.21 @@ -99,9 +99,9 @@
1.22 ; Units for these variables are:
1.23 ; dataob : MgC/ha
1.24 ; datamod0 : gC/m2
1.25 -; We want to convert these to KgC/m2
1.26 +; We want to convert these to kgC/m2
1.27 ; ha = 100m*100m = 10,000 m2
1.28 -; MgC/ha*1000/10,000 = KgC/m2
1.29 +; MgC/ha*1000/10,000 = kgC/m2
1.30
1.31 factor_aboveground = 0.5
1.32 factor_unit_ob = 0.1
1.33 @@ -110,8 +110,8 @@
1.34 dataob = dataob * factor_aboveground * factor_unit_ob
1.35 datamod0 = datamod0 * factor_unit_mod
1.36
1.37 - dataob@units = "KgC/m2"
1.38 - datamod0@units = "KgC/m2"
1.39 + dataob@units = "kg C m~S~-2~N~"
1.40 + datamod0@units = "kg C m~S~-2~N~"
1.41
1.42 dataob@long_name = "Amazon Biomass"
1.43 datamod0@long_name = "Amazon Biomass"
1.44 @@ -139,7 +139,7 @@
1.45 ;********************************************************
1.46 ; sum over amazom_mask area:
1.47
1.48 -; Peta g = 1.e15 g = 1.e12 Kg
1.49 +; Peta g = 1.e15 g = 1.e12 kg
1.50 factor_unit = 1.e-12
1.51
1.52 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
1.53 @@ -204,6 +204,7 @@
1.54
1.55 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.56 "rm "+plot_name+"."+plot_type)
1.57 + ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.58
1.59 ;------------------------------------------------------------------------
1.60 ; contour ob
1.61 @@ -225,6 +226,7 @@
1.62
1.63 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.64 "rm "+plot_name+"."+plot_type)
1.65 + ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.66
1.67 ;------------------------------------------------------------------------
1.68 ; contour model
1.69 @@ -234,7 +236,7 @@
1.70 resg@cnLevelSpacingF = 2. ; interval
1.71
1.72 plot_name = "global_model"
1.73 - title = "Model "+ model_name
1.74 + title = "Model "+ model_name + " (2000)"
1.75 resg@tiMainString = title
1.76
1.77 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.78 @@ -246,6 +248,7 @@
1.79
1.80 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.81 "rm "+plot_name+"."+plot_type)
1.82 + ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.83
1.84 ;------------------------------------------------------------------------
1.85 ; contour model vs ob
1.86 @@ -273,13 +276,17 @@
1.87
1.88 ccrG = esccr(ug,vg,0)
1.89
1.90 - score_max = 5.0
1.91 + score_max = 10.0
1.92
1.93 ; Miomass = (ccrG*ccrG)* score_max
1.94 ; new eq
1.95 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
1.96 + print ("bias=" + bias + "dimsizes(ug)=" + dimsizes(ug))
1.97 + print ("other bias=" + sum(abs(ug-vg)/(ug+vg)))
1.98 Mbiomass = (1. - (bias/dimsizes(ug)))*score_max
1.99 + print ("Mbiomass=" + Mbiomass)
1.100 M_biomass = sprintf("%.2f", Mbiomass)
1.101 + print ("M_biomass=" + M_biomass)
1.102
1.103 if (isvar("compare")) then
1.104 system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \
1.105 @@ -309,7 +316,7 @@
1.106
1.107 ;(b) model
1.108
1.109 - title = "Model "+ model_name
1.110 + title = "Model "+ model_name + " (2000)"
1.111 resg@tiMainString = title
1.112
1.113 plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)
1.114 @@ -318,7 +325,7 @@
1.115
1.116 zz = datamod
1.117 zz = datamod - dataob
1.118 - title = "Model_"+model_name+" - Observed"
1.119 + title = "Model "+model_name+" (2000) - Observed"
1.120
1.121 resg@cnMinLevelValF = -10. ; Min level
1.122 resg@cnMaxLevelValF = 10. ; Max level
1.123 @@ -340,6 +347,7 @@
1.124
1.125 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.126 "rm "+plot_name+"."+plot_type)
1.127 + ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.128
1.129 resg@gsnFrame = True ; draw plot
1.130 resg@gsnDraw = True ; advance frame
1.131 @@ -363,7 +371,7 @@
1.132 gRes@txFontHeightF = 0.02
1.133 gRes@txAngleF = 0
1.134
1.135 - area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" Kg C/m2)"
1.136 + area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" kg C m~S~-2~N~)"
1.137
1.138 gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
1.139 ;-----------------------------------------
1.140 @@ -377,6 +385,7 @@
1.141
1.142 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.143 "rm "+plot_name+"."+plot_type)
1.144 + ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.145
1.146 ;------------------------------------------------------------------------
1.147 ; contour model: masked
1.148 @@ -386,7 +395,7 @@
1.149 resg@cnLevelSpacingF = 2. ; interval
1.150
1.151 plot_name = "global_mask_model"
1.152 - title = "Model "+ model_name
1.153 + title = "Model "+ model_name + " (2000)"
1.154 resg@tiMainString = title
1.155
1.156 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.157 @@ -398,7 +407,7 @@
1.158 gRes@txFontHeightF = 0.02
1.159 gRes@txAngleF = 0
1.160
1.161 - area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" Kg C/m2)"
1.162 + area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" kg C m~S~-2~N~)"
1.163
1.164 gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
1.165 ;-----------------------------------------
1.166 @@ -412,6 +421,7 @@
1.167
1.168 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.169 "rm "+plot_name+"."+plot_type)
1.170 + ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.171
1.172 ;------------------------------------------------------------------------
1.173 ; contour model vs ob: masked
1.174 @@ -434,7 +444,7 @@
1.175 delete (ug)
1.176 delete (vg)
1.177
1.178 - score_max = 5.
1.179 + score_max = 0.
1.180
1.181 uu = ndtooned(zm)
1.182 vv = ndtooned(zo)
1.183 @@ -488,7 +498,7 @@
1.184
1.185 ;(b) model
1.186
1.187 - title = "Model "+ model_name
1.188 + title = "Model "+ model_name + " (2000)"
1.189 resg@tiMainString = title
1.190
1.191 plot(1) = gsn_csm_contour_map_ce(wks,zm,resg)
1.192 @@ -497,7 +507,7 @@
1.193
1.194 zz = zo
1.195 zz = zm - zo
1.196 - title = "Model_"+model_name+" - Observed"
1.197 + title = "Model "+model_name+" (2000) - Observed"
1.198
1.199 resg@cnMinLevelValF = -10. ; Min level
1.200 resg@cnMaxLevelValF = 10. ; Max level
1.201 @@ -518,6 +528,7 @@
1.202
1.203 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.204 "rm "+plot_name+"."+plot_type)
1.205 + ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.206
1.207 ;***************************************************************************
1.208 ; add total score and write to file