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

 ; interpolate into model grids (T31)

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

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

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

 begin

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

 ; read in observed data

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

   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"

   fili  = "amazon_biomass.nc"

   g     = addfile (diri+fili,"r")

   bi    = g->BIOMASS   

   xi    = g->lon 

   yi    = g->lat

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

 ; read in model data

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

  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

  fili2  = "b30.061n_1995-2004_ANN_climo_lnd.nc"

  f      = addfile (diri2+fili2,"r")

  lon    = f->lon     

  lat    = f->lat

  nlon   = dimsizes(lon)

  nlat   = dimsizes(lat)      

 ; print (xi)

 ; print (yi)

  ind_lonL = min(ind(lon .ge. xi(0)))

  ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1)))

 print (ind_lonL)

 print (xi(0))

 print (lon(ind_lonL))

 print (ind_lonR)

 print (xi(dimsizes(xi)-1))

 print (lon(ind_lonR))

  ind_latS = min(ind(lat .ge. yi(0)))

  ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1)))

 print (ind_latS)

 print (yi(0))

 print (lat(ind_latS))

 print (ind_latN)

 print (yi(dimsizes(yi)-1))

 print (lat(ind_latN))

 exit

  bo = new((/nlat,nlon/),float)

     do j=0,nlat-1

        if (j.eq.0 .or. j.eq.nlat-1) then

           if (j.eq.0) then

              LATS = -90.          

              LATN = lat(j)+0.5*(lat(j+1)-lat(j))

           end if

           if (j.eq.nlat-1) then

              LATS = lat(j)-0.5*(lat(j)-lat(j-1))

              LATN = 90.                  

           end if

        else

           LATS = lat(j)-0.5*(lat(j)-lat(j-1))

           LATN = lat(j)+0.5*(lat(j+1)-lat(j))

        end if

 ;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency

 ;      TEMP = bi({LATS:LATN},:)      ; 2D [lat,lon]

       do i=0,nlon-1

        if (i.eq.0 .or. i.eq.nlon-1) then

           if (i.eq.0) then

              LONL = 0.          

              LONR = lon(i)+0.5*(lon(i+1)-lon(i))

           end if

           if (i.eq.nlon-1) then

              LONL = lon(i)-0.5*(lon(i)-lon(i-1))

              LONR = 360.                 

           end if

        else

           LONL = lon(i)-0.5*(lon(i)-lon(i-1))

           LONR = lon(i)+0.5*(lon(i+1)-lon(i))

        end if

 ;print (LATS)

 ;print (LATN)

 ;print (LONL)

 ;print (LONR)

          bo(j,i) = avg(bi({LATS:LATN},{LONL:LONR}))  

 ;        bo(j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)

       end do

 ;     delete(CLAT)

 ;     delete(TEMP) 

     end do

   bo!0   = "lat"

   bo!1   = "lon"

   bo&lat =  lat

   bo&lon =  lon

   bo@units      = bi@units

   bo@long_name  = bi@long_name

 ; bo@_FillValue = bi@_FillValue

   bo@_FillValue = 1.e+36

   c->NPP  = bo

 end