|
1 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
4 |
|
5 begin |
|
6 ;---------------------------------------------------------- |
|
7 year = "2004" |
|
8 |
|
9 nlat = 3600 |
|
10 mlon = 7200 |
|
11 month_in_year = 12 |
|
12 ;---------------------------------------------------------- |
|
13 diri = "/fis/cgd/cseg/people/jeff/clamp_data/" |
|
14 diro = "/fis/cgd/cseg/people/jeff/clamp_data/" |
|
15 filo = "MOD15A2_LAI_" + year + "_monthly.nc" |
|
16 c = addfile(diro+filo,"c") |
|
17 filedimdef(c,"time",-1,True) |
|
18 |
|
19 day = (/"001","009","017","025","033","041","049","057","065","073", \ |
|
20 "081","089","097","105","113","121","129","137","145","153", \ |
|
21 "161","169","177","185","193","201","209","217","225","233", \ |
|
22 "241","249","257","265","273","281","289","297","305","313", \ |
|
23 "321","329","337","345","353","361"/) |
|
24 |
|
25 file_index_L = (/0,3,7 ,11,15,18,22,26,30,34,38,41/) |
|
26 file_index_R = (/3,7,11,14,18,22,26,30,34,37,41,45/) |
|
27 scale_L = (/8,1,5,6,8,1,3,4,5,7,8,2/) |
|
28 scale_R = (/7,3,2,8,7,5,4,3,1,8,6,5/) |
|
29 |
|
30 ; dayI = stringtointeger(day) |
|
31 ; print (dayI) |
|
32 |
|
33 x = new((/12,nlat,mlon/),float) |
|
34 x = 0. |
|
35 x@_FillValue = 1.e+36 |
|
36 |
|
37 do m = 0,month_in_year-1 |
|
38 nday = 0 ; number of day in a month |
|
39 do n = file_index_L(m),file_index_R(m) |
|
40 fili = "MOD15A2_GEO_" + year + day(n) + ".hdf" |
|
41 print (fili) |
|
42 a = addfile(diri+fili,"r") |
|
43 y = a->Lai_0_05deg |
|
44 |
|
45 y@_FillValue = inttobyte(253) ; special missing value? |
|
46 y@_FillValue = inttobyte(255) ; missing value |
|
47 |
|
48 z = byte2flt(y) |
|
49 |
|
50 z@_FillValue = 1.e+36 |
|
51 |
|
52 ; weighted by day, data is 8-day average |
|
53 scale = 8 |
|
54 if (n .eq. file_index_L(m)) then |
|
55 scale = scale_L(m) |
|
56 end if |
|
57 if (n .eq. file_index_R(m)) then |
|
58 scale = scale_R(m) |
|
59 end if |
|
60 |
|
61 nday = nday + scale |
|
62 |
|
63 x(m,:,:) = where(ismissing(z),1.e+36,x(m,:,:)+z(:,:)*scale) |
|
64 |
|
65 delete (a) |
|
66 delete (y) |
|
67 delete (z) |
|
68 end do |
|
69 x(m,:,:) = where(ismissing(x(m,:,:)),1.e+36,x(m,:,:)/nday) |
|
70 print (nday) |
|
71 print (min(x) + "/" + max(x)) |
|
72 end do |
|
73 |
|
74 lat = latGlobeFo(nlat, "lat", "latitude", "degrees_north") |
|
75 lat = (/ lat(::-1) /) ; make N->S |
|
76 lon = lonGlobeFo(mlon, "lon", "longitude", "degrees_east") |
|
77 lon = (/ lon - 180. /) ; subtract 180 from all values |
|
78 lon&lon = lon ; update coordinates |
|
79 |
|
80 x!0 = "time" |
|
81 x!1 = "lat" |
|
82 x!2 = "lon" |
|
83 x&time= ispan(1,month_in_year,1) |
|
84 x&lat= lat |
|
85 x&lon= lon |
|
86 x@units = "none" |
|
87 x@long_name = "Leaf Area Index" |
|
88 |
|
89 c->LAI = x |
|
90 end |
|
91 |