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