diff -r 000000000000 -r 0c6405ab2ff4 biomass/31_area_sum.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/biomass/31_area_sum.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,194 @@
+;*************************************************
+; ce_1.ncl
+;************************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
+;************************************************
+begin
+
+ plot_type     = "ps"
+ plot_type_new = "png"
+
+;***********************************************
+
+ model = "cn"
+;model = "casa"
+
+ model_grid = "T42"
+
+ model_name = "10" + model
+
+;************************************************
+; read amazon mask data
+;************************************************
+
+  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
+  fili  = "amazon_mask_"+ model_grid + ".nc"
+  f     = addfile (diri+fili,"r")
+
+  mask_amazon0 = f->mask_amazon
+
+  delete (f)
+
+;--------------------------------------------------
+; get model data: landfrac and area
+
+  dirs = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/" 
+  fils = "lnd_"+ model_grid +".nc"
+  fs   = addfile (dirs+fils,"r")
+  
+  landfrac = fs->landfrac
+  area0    = fs->area
+
+  delete (fs)
+
+; change area from km**2 to m**2
+  area0 = area0 * 1.e6
+             
+;-----------------------------
+; take into account landfrac
+
+  area0     = area0 * landfrac
+
+;--------------------------------------------
+; read model data
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = "i01.10"+ model +"_1948-2004_ANN_climo.nc"
+ fm   = addfile (dirm+film,"r")
+
+ if (model .eq. "cn") then
+;   unit: gC/m2  
+    data1  = fm->LIVESTEMC
+    data2  = fm->DEADSTEMC
+    data3  = fm->LEAFC
+
+    datamod0 = data1(0,:,:)
+    datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
+ end if
+
+ if (model .eq. "casa") then
+;   unit: gC/m2
+;   factor_WOODC = 0.7  
+    data1  = fm->WOODC
+    data2  = fm->LEAFC
+
+    datamod0 = data1(0,:,:)
+    datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
+ end if
+;----------------------------------------------------
+
+ xm       = fm->lon  
+ ym       = fm->lat
+
+ delete (fm)
+;************************************************
+; read observed data
+;************************************************
+ ob_name = "LC15_Amazon_Biomass"
+
+ diro = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
+ filo = "amazon_biomass_"+model_grid+".nc"
+ fo   = addfile (diro+filo,"r")
+ 
+ dataob   = fo->BIOMASS
+ xo       = fo->lon  
+ yo       = fo->lat
+
+ delete (fo)
+;************************************************
+; Units for these variables are:
+
+; dataob   : MgC/ha
+; datamod0 : gC/m2
+
+; We want to convert these to KgC/m2
+; ha = 100m*100m = 10,000 m2
+; MgC/ha*1000/10,000 = KgC/m2
+; Peta g = 1.e15 g = 1.e12 Kg
+
+  factor_aboveground = 0.5
+  factor_unit_ob     = 0.1 
+  factor_unit_mod    = 0.001          
+
+  dataob   = dataob * factor_aboveground * factor_unit_ob
+  datamod0 = datamod0 * factor_unit_mod 
+
+  dataob@units       = "KgC/m2"
+  datamod0@units     = "KgC/m2"
+
+  dataob@long_name   = "Amazon Biomass"
+  datamod0@long_name = "Amazon Biomass"
+
+;********************************************************
+; get subset of datamod that match dataob
+  
+  nlon = dimsizes(xo)
+  nlat = dimsizes(yo)
+
+  ind_lonL = ind(xm .eq. xo(0))
+  ind_lonR = ind(xm .eq. xo(nlon-1))
+  ind_latS = ind(ym .eq. yo(0))
+  ind_latN = ind(ym .eq. yo(nlat-1))
+
+  datamod  = dataob
+  datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
+
+  area  = dataob
+  area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)
+
+  mask_amazon  = dataob
+  mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)
+
+;********************************************************
+; sum over amazom_mask area:
+
+; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
+
+  area_sum = sum(area*mask_amazon)
+
+  area_sum_amazon_ob  = sum(dataob*area*mask_amazon)  
+  area_sum_amazon_mod = sum(datamod*area*mask_amazon)
+
+; Peta g = 1.e15 g = 1.e12 Kg
+
+  print (area_sum_amazon_ob)
+  print (area_sum_amazon_mod)
+  print (area_sum)
+exit
+;---------------------------------------------------------------------- 
+; global contour 
+
+;res
+  resg                     = True             ; Use plot options
+  resg@cnFillOn            = True             ; Turn on color fill
+  resg@gsnSpreadColors     = True             ; use full colormap
+; resg@cnFillMode          = "RasterFill"     ; Turn on raster color
+; resg@lbLabelAutoStride   = True
+  resg@cnLinesOn           = False            ; Turn off contourn lines
+  resg@mpFillOn            = False            ; Turn off map fill
+  resg@gsnAddCyclic        = False
+
+  resg@gsnSpreadColors      = True            ; use full colormap
+  resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
+  resg@cnMinLevelValF       = 0.              ; Min level
+  resg@cnMaxLevelValF       = 30.             ; Max level
+  resg@cnLevelSpacingF      = 2.              ; interval
+
+  resg@mpMinLatF            = -21.1      ; range to zoom in on
+  resg@mpMaxLatF            =  13.8
+  resg@mpMinLonF            =  277.28
+  resg@mpMaxLonF            =  326.43
+;------------------------------------------------------------------------
+;global contour ob
+  
+  plot_name = "global_ob"
+  title     = ob_name+" "+ model_grid
+  resg@tiMainString  = title
+
+  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
+  gsn_define_colormap(wks,"gui_default")     ; choose colormap
+
+  plot = gsn_csm_contour_map_ce(wks,dataob,resg)   
+  frame(wks) 
+end
\ No newline at end of file