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 = "2004"
forrest@0:   
forrest@0:   nlat  = 3600
forrest@0:   mlon  = 7200
forrest@0:   month_in_year = 12
forrest@0: ;----------------------------------------------------------
forrest@0:   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/"
forrest@0:   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/"
forrest@0:   filo  = "MOD15A2_LAI_" + year + "_monthly.nc"
forrest@0:   c = addfile(diro+filo,"c")
forrest@0:   filedimdef(c,"time",-1,True)
forrest@0: 
forrest@0:   day = (/"001","009","017","025","033","041","049","057","065","073", \
forrest@0:           "081","089","097","105","113","121","129","137","145","153", \
forrest@0:           "161","169","177","185","193","201","209","217","225","233", \
forrest@0:           "241","249","257","265","273","281","289","297","305","313", \
forrest@0:           "321","329","337","345","353","361"/)
forrest@0: 
forrest@0:   file_index_L = (/0,3,7 ,11,15,18,22,26,30,34,38,41/)
forrest@0:   file_index_R = (/3,7,11,14,18,22,26,30,34,37,41,45/)
forrest@0:   scale_L      = (/8,1,5,6,8,1,3,4,5,7,8,2/)
forrest@0:   scale_R      = (/7,3,2,8,7,5,4,3,1,8,6,5/)
forrest@0: 
forrest@0: ;  dayI = stringtointeger(day)
forrest@0: ;  print (dayI)
forrest@0: 
forrest@0:   x = new((/12,nlat,mlon/),float)
forrest@0:   x = 0.
forrest@0:   x@_FillValue =  1.e+36
forrest@0: 
forrest@0:  do m = 0,month_in_year-1
forrest@0:     nday = 0                                   ; number of day in a month
forrest@0:  do n = file_index_L(m),file_index_R(m)
forrest@0:     fili  = "MOD15A2_GEO_" + year + day(n) + ".hdf"
forrest@0:     print (fili)
forrest@0:     a = addfile(diri+fili,"r")
forrest@0:     y = a->Lai_0_05deg
forrest@0: 
forrest@0:     y@_FillValue = inttobyte(253)              ; special missing value? 
forrest@0:     y@_FillValue = inttobyte(255)              ; missing value 
forrest@0: 
forrest@0:     z = byte2flt(y)
forrest@0:           
forrest@0:     z@_FillValue =  1.e+36
forrest@0: 
forrest@0: ;   weighted by day, data is 8-day average                 
forrest@0:     scale = 8              
forrest@0:     if (n .eq. file_index_L(m)) then
forrest@0:        scale = scale_L(m)
forrest@0:     end if
forrest@0:     if (n .eq. file_index_R(m)) then
forrest@0:        scale = scale_R(m)
forrest@0:     end if
forrest@0: 
forrest@0:     nday = nday + scale
forrest@0: 
forrest@0:     x(m,:,:) = where(ismissing(z),1.e+36,x(m,:,:)+z(:,:)*scale)
forrest@0: 
forrest@0:     delete (a)
forrest@0:     delete (y)
forrest@0:     delete (z)
forrest@0:  end do
forrest@0:     x(m,:,:) = where(ismissing(x(m,:,:)),1.e+36,x(m,:,:)/nday)
forrest@0:     print (nday)
forrest@0:     print (min(x) + "/" + max(x))
forrest@0:  end do
forrest@0: 
forrest@0:   lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
forrest@0:   lat  = (/ lat(::-1) /)                     ; make N->S
forrest@0:   lon  = lonGlobeFo(mlon, "lon", "longitude", "degrees_east")
forrest@0:   lon  = (/ lon - 180. /) ; subtract 180 from all values 
forrest@0:   lon&lon = lon           ; update coordinates
forrest@0: 
forrest@0:   x!0  = "time"
forrest@0:   x!1  = "lat"
forrest@0:   x!2  = "lon"
forrest@0:   x&time= ispan(1,month_in_year,1)
forrest@0:   x&lat=  lat
forrest@0:   x&lon=  lon
forrest@0:   x@units      = "none"
forrest@0:   x@long_name  = "Leaf Area Index"
forrest@0: 
forrest@0:   c->LAI  = x 
forrest@0: end
forrest@0: