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

; 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