|
1 ;******************************************************** |
|
2 ; histogram normalized by rain and compute correleration |
|
3 ;******************************************************** |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test" |
|
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test" |
|
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
7 |
|
8 procedure pminmax(data:numeric,name:string) |
|
9 begin |
|
10 print ("min/max " + name + " = " + min(data) + "/" + max(data)) |
|
11 if(isatt(data,"units")) then |
|
12 print (name + " units = " + data@units) |
|
13 end if |
|
14 end |
|
15 |
|
16 ; Main code. |
|
17 begin |
|
18 |
|
19 nclass = 20 |
|
20 |
|
21 plot_type = "ps" |
|
22 plot_type_new = "png" |
|
23 |
|
24 ;************************************************ |
|
25 ; read data: model |
|
26 ;************************************************ |
|
27 co2_i = 283.1878 |
|
28 co2_f = 364.1252 |
|
29 |
|
30 model_grid = "T42" |
|
31 |
|
32 ;model_name_i = "i01.07cn" |
|
33 ;model_name_f = "i01.10cn" |
|
34 |
|
35 model_name_i = "i01.07casa" |
|
36 model_name_f = "i01.10casa" |
|
37 |
|
38 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
39 film_i = model_name_i + "_1990-2004_ANN_climo.nc" |
|
40 film_f = model_name_f + "_1990-2004_ANN_climo.nc" |
|
41 |
|
42 fm_i = addfile (dirm+film_i,"r") |
|
43 fm_f = addfile (dirm+film_f,"r") |
|
44 |
|
45 npp_i = fm_i->NPP |
|
46 npp_f = fm_f->NPP |
|
47 |
|
48 ;************************************************ |
|
49 ; read data: observed |
|
50 ;************************************************ |
|
51 |
|
52 ob_name = "MODIS MOD 15A2 2000-2005" |
|
53 |
|
54 diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" |
|
55 filo = "land_class_"+model_grid+".nc" |
|
56 |
|
57 fo = addfile(diro+filo,"r") |
|
58 |
|
59 classob = tofloat(fo->LAND_CLASS) |
|
60 |
|
61 class_name = (/"Water Bodies" \ |
|
62 ,"Evergreen Needleleaf Forests" \ |
|
63 ,"Evergreen Broadleaf Forests" \ |
|
64 ,"Deciduous Needleleaf Forest" \ |
|
65 ,"Deciduous Broadleaf Forests" \ |
|
66 ,"Mixed Forests" \ |
|
67 ,"Closed Bushlands" \ |
|
68 ,"Open Bushlands" \ |
|
69 ,"Woody Savannas (S. Hem.)" \ |
|
70 ,"Savannas (S. Hem.)" \ |
|
71 ,"Grasslands" \ |
|
72 ,"Permanent Wetlands" \ |
|
73 ,"Croplands" \ |
|
74 ,"Urban and Built-Up" \ |
|
75 ,"Cropland/Natural Vegetation Mosaic" \ |
|
76 ,"Permanent Snow and Ice" \ |
|
77 ,"Barren or Sparsely Vegetated" \ |
|
78 ,"Unclassified" \ |
|
79 ,"Woody Savannas (N. Hem.)" \ |
|
80 ,"Savannas (N. Hem.)" \ |
|
81 /) |
|
82 |
|
83 ;******************************************************************* |
|
84 ; Calculate "nice" bins for binning the data in equally spaced ranges |
|
85 ;******************************************************************** |
|
86 nclassn = nclass + 1 |
|
87 range = fspan(0,nclassn-1,nclassn) |
|
88 ; print (range) |
|
89 |
|
90 ; Use this range information to grab all the values in a |
|
91 ; particular range, and then take an average. |
|
92 |
|
93 nr = dimsizes(range) |
|
94 nx = nr-1 |
|
95 xvalues = new((/2,nx/),float) |
|
96 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. |
|
97 dx = xvalues(0,1) - xvalues(0,0) ; range width |
|
98 dx4 = dx/4 ; 1/4 of the range |
|
99 xvalues(1,:) = xvalues(0,:) - dx/5. |
|
100 |
|
101 ; get data |
|
102 |
|
103 DATA11_1D = ndtooned(classob) |
|
104 DATA12_1D = ndtooned(npp_i) |
|
105 DATA22_1D = ndtooned(npp_f) |
|
106 |
|
107 yvalues = new((/2,nx/),float) |
|
108 mn_yvalues = new((/2,nx/),float) |
|
109 mx_yvalues = new((/2,nx/),float) |
|
110 |
|
111 do nd=0,1 |
|
112 |
|
113 ; See if we are doing model or observational data. |
|
114 |
|
115 if(nd.eq.0) then |
|
116 data_ob = DATA11_1D |
|
117 data_mod = DATA12_1D |
|
118 else |
|
119 data_ob = DATA11_1D |
|
120 data_mod = DATA22_1D |
|
121 end if |
|
122 |
|
123 ; Loop through each range and check for values. |
|
124 |
|
125 do i=0,nr-2 |
|
126 if (i.ne.(nr-2)) then |
|
127 ; print("") |
|
128 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
129 idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1))) |
|
130 else |
|
131 ; print("") |
|
132 ; print("In range ["+range(i)+",)") |
|
133 idx = ind(data_ob.ge.range(i)) |
|
134 end if |
|
135 |
|
136 ; Calculate average, and get min and max. |
|
137 |
|
138 if(.not.any(ismissing(idx))) then |
|
139 yvalues(nd,i) = avg(data_mod(idx)) |
|
140 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
141 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
142 count = dimsizes(idx) |
|
143 else |
|
144 count = 0 |
|
145 yvalues(nd,i) = yvalues@_FillValue |
|
146 mn_yvalues(nd,i) = yvalues@_FillValue |
|
147 mx_yvalues(nd,i) = yvalues@_FillValue |
|
148 end if |
|
149 |
|
150 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
151 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
152 |
|
153 ; Clean up for next time in loop. |
|
154 |
|
155 delete(idx) |
|
156 end do |
|
157 delete(data_ob) |
|
158 delete(data_mod) |
|
159 end do |
|
160 |
|
161 ;============================ |
|
162 ;compute beta |
|
163 ;============================ |
|
164 |
|
165 u = yvalues(0,:) |
|
166 v = yvalues(1,:) |
|
167 |
|
168 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
169 uu = u(good) |
|
170 vv = v(good) |
|
171 ww = class_name(good) |
|
172 |
|
173 n_biome = dimsizes(uu) |
|
174 |
|
175 beta_biome = new((/n_biome/),float) |
|
176 |
|
177 beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i) |
|
178 |
|
179 beta_biome_avg = avg(beta_biome) |
|
180 |
|
181 print("class/beta: " + ww + "/" + beta_biome) |
|
182 print (beta_biome_avg) |
|
183 |
|
184 end |
|
185 |