forrest@0: ; ***********************************************
forrest@0: ; interpolate into model grids (T31)
forrest@0: ; ***********************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;************************************************
forrest@0: begin
forrest@0: 
forrest@0:   year = 2003
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; output data
forrest@0: ;************************************************
forrest@0:   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
forrest@0:   filo  = "LAI_" + year + "_monthly_1.9.nc"
forrest@0: 
forrest@0:   c = addfile(diro+filo,"c")
forrest@0:   filedimdef(c,"time",-1,True)
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; read in observed data
forrest@0: ;************************************************
forrest@0:   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
forrest@0:   fili  = "LAI_" + year + "_monthly.nc"
forrest@0: 
forrest@0:   g     = addfile (diri+fili,"r")
forrest@0:   bi    = g->LAI   
forrest@0:   xi    = g->lon 
forrest@0:   yi    = g->lat
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; change into 0-360E, 90S-90N
forrest@0: ;************************************************
forrest@0:  
forrest@0:  yi    = (/ yi(::-1) /)
forrest@0:  bi    = (/ bi(:,::-1,:) /)
forrest@0:  printVarSummary(bi)
forrest@0: 
forrest@0:  b2    = bi
forrest@0:  x2    = xi   
forrest@0:  
forrest@0:  nx = dimsizes(xi)
forrest@0:  do i= 0,nx-1
forrest@0:     if (i .lt. 3600) then
forrest@0:        p = i + 3600
forrest@0:        xi(p) = x2(i) + 360.      
forrest@0:     else
forrest@0:        p = i - 3600
forrest@0:        xi(p) = x2(i)
forrest@0:     end if
forrest@0:     bi(:,:,p)= b2(:,:,i) 
forrest@0:  end do
forrest@0: 
forrest@0:  bi&lat =  yi
forrest@0:  bi&lon =  xi
forrest@0: 
forrest@0: ;print (xi)
forrest@0: ;print (yi)
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; read in model data
forrest@0: ;************************************************
forrest@0:  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0:  fili2  = "newcn05_ncep_1i_MONS_climo_lnd.nc"
forrest@0:  f      = addfile (diri2+fili2,"r")
forrest@0: 
forrest@0:  lon    = f->lon     
forrest@0:  lat    = f->lat
forrest@0:  nlon   = dimsizes(lon)
forrest@0:  nlat   = dimsizes(lat)      
forrest@0: 
forrest@0: ; print (xi)
forrest@0: ; print (yi)
forrest@0: ; print (xo)
forrest@0: ; print (yo)
forrest@0: 
forrest@0:  bo = new((/12,nlat,nlon/),float)
forrest@0: 
forrest@0:  do m = 0,11
forrest@0:     do j=0,nlat-1
forrest@0:        if (j.eq.0 .or. j.eq.nlat-1) then
forrest@0:           if (j.eq.0) then
forrest@0:              LATS = -90.          
forrest@0:              LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0:           end if
forrest@0:           if (j.eq.nlat-1) then
forrest@0:              LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0:              LATN = 90.                  
forrest@0:           end if
forrest@0:        else
forrest@0:           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
forrest@0:           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
forrest@0:        end if
forrest@0:  
forrest@0:       do i=0,nlon-1
forrest@0:        if (i.eq.0 .or. i.eq.nlon-1) then
forrest@0:           if (i.eq.0) then
forrest@0:              LONL = 0.          
forrest@0:              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0:           end if
forrest@0:           if (i.eq.nlon-1) then
forrest@0:              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0:              LONR = 360.                 
forrest@0:           end if
forrest@0:        else
forrest@0:           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
forrest@0:           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
forrest@0:        end if
forrest@0: 
forrest@0: ;print (LATS)
forrest@0: ;print (LATN)
forrest@0: ;print (LONL)
forrest@0: ;print (LONR)
forrest@0: 
forrest@0:          bo(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))  
forrest@0:       end do 
forrest@0:     end do
forrest@0:   end do
forrest@0: 
forrest@0:   bo!0   = "time"
forrest@0:   bo!1   = "lat"
forrest@0:   bo!2   = "lon"
forrest@0:   bo&time= bi&time
forrest@0:   bo&lat = lat
forrest@0:   bo&lon = lon
forrest@0:   bo@units      = bi@units
forrest@0:   bo@long_name  = bi@long_name
forrest@0:   bo@_FillValue = bi@_FillValue
forrest@0: 
forrest@0:   c->LAI  = bo
forrest@0: end