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

 ; 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

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

 ; output data

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

   diro  = ""

   filo  = "npp_T42_mean_2000-2004.nc"

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

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

 ; read in observed data

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

   diri  = ""

   fili  = "npp_0.05deg_mean_2000-2004.nc"

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

   bi    = g->NPP   

   xi    = g->lon 

   yi    = g->lat

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

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

 ; Observed NPP*scale_factor

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

  ;scale_factor = 0.1

  ; FMH : Now the scaling occurs in the previous routine the writes out the mean

  ;       netCDF file, so don't scale again here.

  scale_factor = 1.

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

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

  b2    = bi

  x2    = xi   

  nx = dimsizes(xi)

  do i= 0,nx-1

     if (i .lt. 3600) then

        p = i + 3600

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

     else

        p = i - 3600

        xi(p) = x2(i)

     end if

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

  end do

  bi&lat =  yi

  bi&lon =  xi

 ;print (xi)

 ;print (yi)

 ;exit

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

 ; read in model data

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

 ;fili2  = "i01.04casa_1605-1629_ANN_climo.nc"

  diri   = "../../model/"

  fili2  = "i01.10cn_2000-2004_ANN_climo.nc"

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

  lon    = f->lon     

  lat    = f->lat

  nlon   = dimsizes(lon)

  nlat   = dimsizes(lat)      

 ; print (xi)

 ; print (yi)

 ; print (xo)

 ; print (yo)

 ;(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!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