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

 ; 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 model grid

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

  model_grid = "T31"

 ;model_grid = "T42"

 ;model_grid = "1.9"

  diri   = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/"

  fili   = "lnd_"+ model_grid +".nc"

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

  lon    = f->lon     

  lat    = f->lat

  nlon   = dimsizes(lon)

  nlat   = dimsizes(lat)      

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

 ; output data

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

   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"

   filo  = "amazon_mask_"+model_grid + ".nc"

   c = addfile(diro+filo,"c")

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

 ; read 1x1 amazon_mask data

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

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

   fili  = "amazon_mask_1x1.nc"

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

   bi    = g->mask_amazon   

   xi    = g->lon 

   yi    = g->lat

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

 ; change from -180-180, 90N-90S into 0-360E, 90S-90N

 ; Observed NPP*scale_factor

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

  yi    = (/ yi(::-1) /)

  bi    =  (/ bi(::-1,:) /)

  b2    = bi

  x2    = xi   

  nx = dimsizes(xi)

  do i= 0,nx-1

     if (i .lt. 180) then

        p = i + 180

        xi(p) = x2(i) + 360.      

     else

        p = i - 180

        xi(p) = x2(i)

     end if

     bi(:,p)= b2(:,i) 

  end do

  bi&lat =  yi

  bi&lon =  xi

 ;print (xi)

 ;print (yi)

 ;(A)

 ;bo = linint2_Wrap(xi,yi,bi,True,xo,yo,0)

 ;(B)

   rad   = 4.*atan(1.)/180.

   clat  = lat

   clat  = cos(lat*rad)

   clat@long_name = "cos(latitude)"

   delete(clat@units)

 ; printVarSummary(clat)

  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 = where(bo.ge.0.5, 1.0,0.)

   bo!0   = "lat"

   bo!1   = "lon"

   bo&lat =  lat

   bo&lon =  lon

 ; bo@units      = bi@units

   bo@long_name  = "amazon mask"

 ; bo@_FillValue = bi@_FillValue

   bo@_FillValue = 1.e+36

   c->mask_amazon  = bo

 end