forrest@2: ; ***********************************************
forrest@2: ; interpolate into model grids (T31)
forrest@2: ; ***********************************************
forrest@2: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@2: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@2: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@2: ;************************************************
forrest@2: begin
forrest@2: ;************************************************
forrest@2: ; output data
forrest@2: ;************************************************
forrest@2:   diro  = ""
forrest@2:   filo  = "npp_T42_mean_2000-2004.nc"
forrest@2:   c = addfile(diro+filo,"c")
forrest@2: 
forrest@2: ;************************************************
forrest@2: ; read in observed data
forrest@2: ;************************************************
forrest@2:   diri  = ""
forrest@2:   fili  = "npp_0.05deg_mean_2000-2004.nc"
forrest@2:   g     = addfile (diri+fili,"r")
forrest@2:   bi    = g->NPP   
forrest@2:   xi    = g->lon 
forrest@2:   yi    = g->lat
forrest@2: 
forrest@2: ;************************************************
forrest@2: ; change into 0-360E, 90S-90N
forrest@2: ; Observed NPP*scale_factor
forrest@2: ;************************************************
forrest@2: 
forrest@2:  ;scale_factor = 0.1
forrest@2:  ; FMH : Now the scaling occurs in the previous routine the writes out the mean
forrest@2:  ;       netCDF file, so don't scale again here.
forrest@2:  scale_factor = 1.
forrest@2:  
forrest@2:  yi    = (/ yi(::-1) /)
forrest@2:  bi    =  (/ bi(::-1,:) /)
forrest@2:  b2    = bi
forrest@2:  x2    = xi   
forrest@2:  
forrest@2:  nx = dimsizes(xi)
forrest@2:  do i= 0,nx-1
forrest@2:     if (i .lt. 3600) then
forrest@2:        p = i + 3600
forrest@2:        xi(p) = x2(i) + 360.      
forrest@2:     else
forrest@2:        p = i - 3600
forrest@2:        xi(p) = x2(i)
forrest@2:     end if
forrest@2:     bi(:,p)= b2(:,i) * scale_factor
forrest@2:  end do
forrest@2: 
forrest@2:  bi&lat =  yi
forrest@2:  bi&lon =  xi
forrest@2: 
forrest@2: ;print (xi)
forrest@2: ;print (yi)
forrest@2: ;exit
forrest@2: 
forrest@2: ;************************************************
forrest@2: ; read in model data
forrest@2: ;************************************************
forrest@2: ;fili2  = "i01.04casa_1605-1629_ANN_climo.nc"
forrest@2:  diri   = "../../model/"
forrest@2:  fili2  = "i01.10cn_2000-2004_ANN_climo.nc"
forrest@2:  f     = addfile (diri+fili2,"r")
forrest@2: 
forrest@2:  lon    = f->lon     
forrest@2:  lat    = f->lat
forrest@2:  nlon   = dimsizes(lon)
forrest@2:  nlat   = dimsizes(lat)      
forrest@2: 
forrest@2: ; print (xi)
forrest@2: ; print (yi)
forrest@2: ; print (xo)
forrest@2: ; print (yo)
forrest@2: 
forrest@2: ;(A)
forrest@2: ;bo = linint2_Wrap(xi,yi,bi,True,xo,yo,0)
forrest@2: 
forrest@2: ;(B)
forrest@2: 
forrest@2:   rad   = 4.*atan(1.)/180.
forrest@2:   clat  = lat
forrest@2:   clat  = cos(lat*rad)
forrest@2:   clat@long_name = "cos(latitude)"
forrest@2:   delete(clat@units)
forrest@2:   printVarSummary(clat)
forrest@2: 
forrest@2:  bo = new((/nlat,nlon/),float)
forrest@2:     do j=0,nlat-1
forrest@2:        if (j.eq.0 .or. j.eq.nlat-1) then
forrest@2:           if (j.eq.0) then
forrest@2:              LATS = -90.          
forrest@2:              LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@2:           end if
forrest@2:           if (j.eq.nlat-1) then
forrest@2:              LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@2:              LATN = 90.                  
forrest@2:           end if
forrest@2:        else
forrest@2:           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@2:           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@2:        end if
forrest@2:  
forrest@2: ;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
forrest@2: ;      TEMP = bi({LATS:LATN},:)      ; 2D [lat,lon]
forrest@2:  
forrest@2:       do i=0,nlon-1
forrest@2:        if (i.eq.0 .or. i.eq.nlon-1) then
forrest@2:           if (i.eq.0) then
forrest@2:              LONL = 0.          
forrest@2:              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@2:           end if
forrest@2:           if (i.eq.nlon-1) then
forrest@2:              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@2:              LONR = 360.                 
forrest@2:           end if
forrest@2:        else
forrest@2:           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@2:           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@2:        end if
forrest@2: 
forrest@2: ;print (LATS)
forrest@2: ;print (LATN)
forrest@2: ;print (LONL)
forrest@2: ;print (LONR)
forrest@2:          bo(j,i) = avg(bi({LATS:LATN},{LONL:LONR}))  
forrest@2: ;        bo(j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
forrest@2:       end do
forrest@2:  
forrest@2: ;     delete(CLAT)
forrest@2: ;     delete(TEMP) 
forrest@2:     end do
forrest@2: 
forrest@2:   bo!0   = "lat"
forrest@2:   bo!1   = "lon"
forrest@2:   bo&lat =  lat
forrest@2:   bo&lon =  lon
forrest@2:   bo@units      = bi@units
forrest@2:   bo@long_name  = bi@long_name
forrest@2: ; bo@_FillValue = bi@_FillValue
forrest@2:   bo@_FillValue = 1.e+36
forrest@2: 
forrest@2:   c->NPP  = bo
forrest@2: end