1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lai/02.read_lai_avg.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,91 @@
     1.4 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
     1.5 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
     1.6 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
     1.7 +
     1.8 +begin 
     1.9 +;----------------------------------------------------------
    1.10 +  year = "2004"
    1.11 +  
    1.12 +  nlat  = 3600
    1.13 +  mlon  = 7200
    1.14 +  month_in_year = 12
    1.15 +;----------------------------------------------------------
    1.16 +  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/"
    1.17 +  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/"
    1.18 +  filo  = "MOD15A2_LAI_" + year + "_monthly.nc"
    1.19 +  c = addfile(diro+filo,"c")
    1.20 +  filedimdef(c,"time",-1,True)
    1.21 +
    1.22 +  day = (/"001","009","017","025","033","041","049","057","065","073", \
    1.23 +          "081","089","097","105","113","121","129","137","145","153", \
    1.24 +          "161","169","177","185","193","201","209","217","225","233", \
    1.25 +          "241","249","257","265","273","281","289","297","305","313", \
    1.26 +          "321","329","337","345","353","361"/)
    1.27 +
    1.28 +  file_index_L = (/0,3,7 ,11,15,18,22,26,30,34,38,41/)
    1.29 +  file_index_R = (/3,7,11,14,18,22,26,30,34,37,41,45/)
    1.30 +  scale_L      = (/8,1,5,6,8,1,3,4,5,7,8,2/)
    1.31 +  scale_R      = (/7,3,2,8,7,5,4,3,1,8,6,5/)
    1.32 +
    1.33 +;  dayI = stringtointeger(day)
    1.34 +;  print (dayI)
    1.35 +
    1.36 +  x = new((/12,nlat,mlon/),float)
    1.37 +  x = 0.
    1.38 +  x@_FillValue =  1.e+36
    1.39 +
    1.40 + do m = 0,month_in_year-1
    1.41 +    nday = 0                                   ; number of day in a month
    1.42 + do n = file_index_L(m),file_index_R(m)
    1.43 +    fili  = "MOD15A2_GEO_" + year + day(n) + ".hdf"
    1.44 +    print (fili)
    1.45 +    a = addfile(diri+fili,"r")
    1.46 +    y = a->Lai_0_05deg
    1.47 +
    1.48 +    y@_FillValue = inttobyte(253)              ; special missing value? 
    1.49 +    y@_FillValue = inttobyte(255)              ; missing value 
    1.50 +
    1.51 +    z = byte2flt(y)
    1.52 +          
    1.53 +    z@_FillValue =  1.e+36
    1.54 +
    1.55 +;   weighted by day, data is 8-day average                 
    1.56 +    scale = 8              
    1.57 +    if (n .eq. file_index_L(m)) then
    1.58 +       scale = scale_L(m)
    1.59 +    end if
    1.60 +    if (n .eq. file_index_R(m)) then
    1.61 +       scale = scale_R(m)
    1.62 +    end if
    1.63 +
    1.64 +    nday = nday + scale
    1.65 +
    1.66 +    x(m,:,:) = where(ismissing(z),1.e+36,x(m,:,:)+z(:,:)*scale)
    1.67 +
    1.68 +    delete (a)
    1.69 +    delete (y)
    1.70 +    delete (z)
    1.71 + end do
    1.72 +    x(m,:,:) = where(ismissing(x(m,:,:)),1.e+36,x(m,:,:)/nday)
    1.73 +    print (nday)
    1.74 +    print (min(x) + "/" + max(x))
    1.75 + end do
    1.76 +
    1.77 +  lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
    1.78 +  lat  = (/ lat(::-1) /)                     ; make N->S
    1.79 +  lon  = lonGlobeFo(mlon, "lon", "longitude", "degrees_east")
    1.80 +  lon  = (/ lon - 180. /) ; subtract 180 from all values 
    1.81 +  lon&lon = lon           ; update coordinates
    1.82 +
    1.83 +  x!0  = "time"
    1.84 +  x!1  = "lat"
    1.85 +  x!2  = "lon"
    1.86 +  x&time= ispan(1,month_in_year,1)
    1.87 +  x&lat=  lat
    1.88 +  x&lon=  lon
    1.89 +  x@units      = "none"
    1.90 +  x@long_name  = "Leaf Area Index"
    1.91 +
    1.92 +  c->LAI  = x 
    1.93 +end
    1.94 +