begin

  filename = "LAI_2000-2005_ANN_0.9.nc"
  
  in_file1 = addfile("surfdata_0.9x1.25_c070501.nc","r")
  in_file2 = addfile("lnd_0.9.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_PFT
  lsmpftdatamask = in_file1->PFTDATA_MASK
  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 latbigindex = 0,latbigdim-1
      do lonbigindex = 0,lonbigdim-1
          if (.not. ismissing(cleanannlai(latbigindex,lonbigindex))) then
          newannlai(latbigindex*2,lonbigindex*2) = cleanannlai(latbigindex,lonbigindex)
          newannlai(latbigindex*2,lonbigindex*2+1) = cleanannlai(latbigindex,lonbigindex)
          newannlai(latbigindex*2+1,lonbigindex*2) = cleanannlai(latbigindex,lonbigindex)
          newannlai(latbigindex*2+1,lonbigindex*2+1) = 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