diff -r 000000000000 -r 0c6405ab2ff4 lai/03.lai_ob_to_0.25deg.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lai/03.lai_ob_to_0.25deg.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,154 @@
+; ***********************************************
+; 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
+
+  year = 2000
+
+;************************************************
+; output data
+;************************************************
+  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
+; filo  = "LAI_" + year + "_monthly_0.25deg.nc"
+  filo  = "LAI_2000-2005_mean_0.25deg.nc"
+; filo  = "LAI_2000-2005_ensemble_0.25deg.nc"
+  c = addfile(diro+filo,"c")
+  filedimdef(c,"time",-1,True)
+
+;************************************************
+; read in observed data
+;************************************************
+  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
+; fili  = "LAI_" + year + "_monthly.nc"
+  fili  = "LAI_2000-2005_mean_0.05deg.nc"
+; fili  = "LAI_2000-2005_ensemble_0.05deg.nc"
+  g     = addfile (diri+fili,"r")
+  bi    = g->LAI   
+  xi    = g->lon 
+  yi    = g->lat
+
+;************************************************
+; change into 0-360E, 90S-90N
+; Observed NPP*scale_factor
+;************************************************
+ 
+ yi    = (/ yi(::-1) /)
+ bi    = (/ bi(:,::-1,:) /)
+ printVarSummary(bi)
+
+ 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) 
+ end do
+
+ bi&lat =  yi
+ bi&lon =  xi
+
+;print (xi)
+;print (yi)
+;exit
+
+;************************************************
+; create 0.25deg lat and lon
+;************************************************
+  nlon   = 360*4
+  nlat   = 180*4
+
+  lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north") ; S->N
+; lat  = (/ lat(::-1) /)                                      ; N->S  
+  lat&lat = lat
+  lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east") ; 0->360 
+; lon  = (/ lon - 180. /) ; subtract 180 from all values      ; 180W-180E
+  lon&lon = lon           ; update coordinates
+
+; print (lon)
+; print (lat)
+
+; rad   = 4.*atan(1.)/180.
+; clat  = lat
+; clat  = cos(lat*rad)
+; clat@long_name = "cos(latitude)"
+; delete(clat@units)
+; printVarSummary(clat)
+
+;bo = new((/12,nlat,nlon/),float)
+ bo = new((/1,nlat,nlon/),float)
+
+;do m = 0,11
+ do m = 0,0
+    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(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))  
+;        bo(m,j,i) = wgt_areaave(TEMP(m,{LONL:LONR}), CLAT, 1.0, 0)
+      end do
+ 
+;     delete(CLAT)
+;     delete(TEMP) 
+    end do
+  end do
+
+  bo!0   = "time"
+  bo!1   = "lat"
+  bo!2   = "lon"
+  bo&time= bi&time
+; bo&time= 1
+  bo&lat = lat
+  bo&lon = lon
+; bo@units      = bi@units
+; bo@long_name  = bi@long_name
+  bo@units      = "none"
+  bo@long_name  = "Leaf Area Index"
+  bo@_FillValue = bi@_FillValue
+
+  c->LAI  = bo
+end
\ No newline at end of file