begin

  filename = "LAI_2000-2005_ANN_0.5.nc"
  
  in_file1 = addfile("surfdata_360x720_mcrop2000.nc","r")
  in_file2 = addfile("lnd_0.5.nc","r")
  in_file3 = addfile("LAI_2000-2005_ANN_1.9.nc","r")
  out_file = addfile(filename,"c")

  fileAttr = True
  fileAttr@title = filename
  fileAttr@creation_date = systemfunc("date")
  fileAttr@source = "Peter Lawrence, NCAR"
  fileattdef(out_file,fileAttr)
  
  lsmedgen = in_file1->EDGEN
  lsmedgee = in_file1->EDGEE
  lsmedges = in_file1->EDGES
  lsmedgew = in_file1->EDGEW
  lsmlatn = in_file1->LATN
  lsmlone = in_file1->LONE
  lsmlats = in_file1->LATS
  lsmlonw = in_file1->LONW
  lsmarea = in_file1->AREA
  lsmnumlon = in_file1->NUMLON
  lsmlongxy = in_file1->LONGXY
  lsmlatixy = in_file1->LATIXY
  lsmlandfracpft = in_file1->LANDFRAC
  lsmpftdatamask = in_file1->LANDMASK
  lsmpft = in_file1->PCT_PFT
  
  lndlandmask = in_file2->landmask
  annlai = in_file3->LAI
    
  latlondims = dimsizes(lsmlongxy)
  latdim = latlondims(0)
  londim = latlondims(1)
  
  latlonbigdims = dimsizes(annlai)
  latbigdim = latlonbigdims(0)
  lonbigdim = latlonbigdims(1)

  dimNames = (/ "lon",  "lat"      /)
  dimSizes = (/ londim, latdim /)
  dimUnlim = (/ False,  False  /)
  filedimdef( out_file, dimNames, dimSizes, dimUnlim )

  filevardef(out_file, "LAI",typeof(lsmpft), (/ "lat", "lon" /))
  filevardef(out_file, "lat",typeof(lsmpft), (/ "lat" /))
  filevardef(out_file, "lon",typeof(lsmpft), (/ "lon" /))

  cleanannlai = where(annlai .ne. annlai@_FillValue,annlai,0.0)
  cleanannlai@missing_value = cleanannlai@_FillValue
  
  newannlai = new((/ latdim, londim /),float)
  newannlai = (/ 0.0 /)
  
  do latindex = 0,latdim-1
      do lonindex = 0,londim-1
      
          latbigindex = floattointeger(latindex * 0.5 / 1.875)
          lonbigindex = floattointeger(lonindex * 0.5 / 2.5)

          if (.not. ismissing(cleanannlai(latbigindex,lonbigindex))) then
              newannlai(latindex,lonindex) = cleanannlai(latbigindex,lonbigindex)
	  end if

      end do
  end do
  
  out_file->LAI@long_name = "Leaf Area Index"
  out_file->LAI@units = "none"
  out_file->LAI@_FillValue = 1.0e36
  out_file->LAI = (/ where(lndlandmask.eq.1,newannlai,out_file->LAI@_FillValue) /)
  out_file->lat = in_file2->lat
  out_file->lon = in_file2->lon
          
end