1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/biomass/31_area_sum.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,194 @@
1.4 +;*************************************************
1.5 +; ce_1.ncl
1.6 +;************************************************
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.9 +;************************************************
1.10 +begin
1.11 +
1.12 + plot_type = "ps"
1.13 + plot_type_new = "png"
1.14 +
1.15 +;***********************************************
1.16 +
1.17 + model = "cn"
1.18 +;model = "casa"
1.19 +
1.20 + model_grid = "T42"
1.21 +
1.22 + model_name = "10" + model
1.23 +
1.24 +;************************************************
1.25 +; read amazon mask data
1.26 +;************************************************
1.27 +
1.28 + diri = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
1.29 + fili = "amazon_mask_"+ model_grid + ".nc"
1.30 + f = addfile (diri+fili,"r")
1.31 +
1.32 + mask_amazon0 = f->mask_amazon
1.33 +
1.34 + delete (f)
1.35 +
1.36 +;--------------------------------------------------
1.37 +; get model data: landfrac and area
1.38 +
1.39 + dirs = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/"
1.40 + fils = "lnd_"+ model_grid +".nc"
1.41 + fs = addfile (dirs+fils,"r")
1.42 +
1.43 + landfrac = fs->landfrac
1.44 + area0 = fs->area
1.45 +
1.46 + delete (fs)
1.47 +
1.48 +; change area from km**2 to m**2
1.49 + area0 = area0 * 1.e6
1.50 +
1.51 +;-----------------------------
1.52 +; take into account landfrac
1.53 +
1.54 + area0 = area0 * landfrac
1.55 +
1.56 +;--------------------------------------------
1.57 +; read model data
1.58 +
1.59 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.60 + film = "i01.10"+ model +"_1948-2004_ANN_climo.nc"
1.61 + fm = addfile (dirm+film,"r")
1.62 +
1.63 + if (model .eq. "cn") then
1.64 +; unit: gC/m2
1.65 + data1 = fm->LIVESTEMC
1.66 + data2 = fm->DEADSTEMC
1.67 + data3 = fm->LEAFC
1.68 +
1.69 + datamod0 = data1(0,:,:)
1.70 + datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
1.71 + end if
1.72 +
1.73 + if (model .eq. "casa") then
1.74 +; unit: gC/m2
1.75 +; factor_WOODC = 0.7
1.76 + data1 = fm->WOODC
1.77 + data2 = fm->LEAFC
1.78 +
1.79 + datamod0 = data1(0,:,:)
1.80 + datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
1.81 + end if
1.82 +;----------------------------------------------------
1.83 +
1.84 + xm = fm->lon
1.85 + ym = fm->lat
1.86 +
1.87 + delete (fm)
1.88 +;************************************************
1.89 +; read observed data
1.90 +;************************************************
1.91 + ob_name = "LC15_Amazon_Biomass"
1.92 +
1.93 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
1.94 + filo = "amazon_biomass_"+model_grid+".nc"
1.95 + fo = addfile (diro+filo,"r")
1.96 +
1.97 + dataob = fo->BIOMASS
1.98 + xo = fo->lon
1.99 + yo = fo->lat
1.100 +
1.101 + delete (fo)
1.102 +;************************************************
1.103 +; Units for these variables are:
1.104 +
1.105 +; dataob : MgC/ha
1.106 +; datamod0 : gC/m2
1.107 +
1.108 +; We want to convert these to KgC/m2
1.109 +; ha = 100m*100m = 10,000 m2
1.110 +; MgC/ha*1000/10,000 = KgC/m2
1.111 +; Peta g = 1.e15 g = 1.e12 Kg
1.112 +
1.113 + factor_aboveground = 0.5
1.114 + factor_unit_ob = 0.1
1.115 + factor_unit_mod = 0.001
1.116 +
1.117 + dataob = dataob * factor_aboveground * factor_unit_ob
1.118 + datamod0 = datamod0 * factor_unit_mod
1.119 +
1.120 + dataob@units = "KgC/m2"
1.121 + datamod0@units = "KgC/m2"
1.122 +
1.123 + dataob@long_name = "Amazon Biomass"
1.124 + datamod0@long_name = "Amazon Biomass"
1.125 +
1.126 +;********************************************************
1.127 +; get subset of datamod that match dataob
1.128 +
1.129 + nlon = dimsizes(xo)
1.130 + nlat = dimsizes(yo)
1.131 +
1.132 + ind_lonL = ind(xm .eq. xo(0))
1.133 + ind_lonR = ind(xm .eq. xo(nlon-1))
1.134 + ind_latS = ind(ym .eq. yo(0))
1.135 + ind_latN = ind(ym .eq. yo(nlat-1))
1.136 +
1.137 + datamod = dataob
1.138 + datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
1.139 +
1.140 + area = dataob
1.141 + area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)
1.142 +
1.143 + mask_amazon = dataob
1.144 + mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)
1.145 +
1.146 +;********************************************************
1.147 +; sum over amazom_mask area:
1.148 +
1.149 +; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
1.150 +
1.151 + area_sum = sum(area*mask_amazon)
1.152 +
1.153 + area_sum_amazon_ob = sum(dataob*area*mask_amazon)
1.154 + area_sum_amazon_mod = sum(datamod*area*mask_amazon)
1.155 +
1.156 +; Peta g = 1.e15 g = 1.e12 Kg
1.157 +
1.158 + print (area_sum_amazon_ob)
1.159 + print (area_sum_amazon_mod)
1.160 + print (area_sum)
1.161 +exit
1.162 +;----------------------------------------------------------------------
1.163 +; global contour
1.164 +
1.165 +;res
1.166 + resg = True ; Use plot options
1.167 + resg@cnFillOn = True ; Turn on color fill
1.168 + resg@gsnSpreadColors = True ; use full colormap
1.169 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.170 +; resg@lbLabelAutoStride = True
1.171 + resg@cnLinesOn = False ; Turn off contourn lines
1.172 + resg@mpFillOn = False ; Turn off map fill
1.173 + resg@gsnAddCyclic = False
1.174 +
1.175 + resg@gsnSpreadColors = True ; use full colormap
1.176 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.177 + resg@cnMinLevelValF = 0. ; Min level
1.178 + resg@cnMaxLevelValF = 30. ; Max level
1.179 + resg@cnLevelSpacingF = 2. ; interval
1.180 +
1.181 + resg@mpMinLatF = -21.1 ; range to zoom in on
1.182 + resg@mpMaxLatF = 13.8
1.183 + resg@mpMinLonF = 277.28
1.184 + resg@mpMaxLonF = 326.43
1.185 +;------------------------------------------------------------------------
1.186 +;global contour ob
1.187 +
1.188 + plot_name = "global_ob"
1.189 + title = ob_name+" "+ model_grid
1.190 + resg@tiMainString = title
1.191 +
1.192 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.193 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.194 +
1.195 + plot = gsn_csm_contour_map_ce(wks,dataob,resg)
1.196 + frame(wks)
1.197 +end
1.198 \ No newline at end of file