;*************************************************

 ; 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