|
1 |
|
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
4 |
|
5 procedure pminmax(data:numeric,name:string) |
|
6 begin |
|
7 print ("min/max " + name + " = " + min(data) + "/" + max(data)) |
|
8 if(isatt(data,"units")) then |
|
9 print (name + " units = " + data@units) |
|
10 end if |
|
11 end |
|
12 |
|
13 ; |
|
14 ; Main code. |
|
15 ; |
|
16 begin |
|
17 data_types = (/ "Obs", "Model" /) |
|
18 data_names = (/ "data.81.nc", "i01.03cn_1545-1569_ANN_climo.nc" /) |
|
19 ; data_names = (/ "data.81.nc", "i01.04casa_1605-1629_ANN_climo.nc" /) |
|
20 ; filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /) |
|
21 filevar_names = (/ (/"PREC_ANN","ANPP_C"/), (/"RAIN","AGNPP"/) /) |
|
22 ndata_types = dimsizes(data_types) |
|
23 |
|
24 data_file_obs = addfile(data_names(0),"r") ; Open obs file |
|
25 data_file_mod = addfile(data_names(1),"r") ; Open model file |
|
26 |
|
27 ;************************************************ |
|
28 ; read in data: observed |
|
29 ;************************************************ |
|
30 PREC_ANN = tofloat(data_file_obs->PREC_ANN) |
|
31 TNPP_C = data_file_obs->TNPP_C |
|
32 xo = data_file_obs->LONG_DD |
|
33 yo = data_file_obs->LAT_DD |
|
34 |
|
35 ; change longitude from (-180,180) to (0,360) |
|
36 ;print (xo) |
|
37 nx = dimsizes(xo) |
|
38 do i= 0,nx-1 |
|
39 if (xo(i) .lt. 0.) then |
|
40 xo(i) = xo(i)+ 360. |
|
41 end if |
|
42 end do |
|
43 ;print (xo) |
|
44 |
|
45 ;************************************************ |
|
46 ; read in data: model |
|
47 ;************************************************ |
|
48 ai = data_file_mod->RAIN |
|
49 bi = data_file_mod->NPP |
|
50 xi = data_file_mod->lon |
|
51 yi = data_file_mod->lat |
|
52 |
|
53 ;************************************************ |
|
54 ; interpolate from model grid to observed grid |
|
55 ;************************************************ |
|
56 RAIN = linint2_points(xi,yi,ai,True,xo,yo,0) |
|
57 NPP = linint2_points(xi,yi,bi,True,xo,yo,0) |
|
58 |
|
59 ;************************************************ |
|
60 ; convert unit |
|
61 ;************************************************ |
|
62 ; Units for these four variables are: |
|
63 ; |
|
64 ; PREC_ANN : mm/year |
|
65 ; TNPP_C : g C/m^2/year |
|
66 ; RAIN : mm/s |
|
67 ; NPP : g C/m^2/s |
|
68 ; |
|
69 ; We want to convert these to "m/year" and "g C/m^2/year". |
|
70 ; |
|
71 nsec_per_year = 60*60*24*365 ; # seconds per year |
|
72 |
|
73 ; Do the necessary conversions. |
|
74 PREC_ANN = PREC_ANN / 1000. |
|
75 RAIN = (RAIN / 1000.) * nsec_per_year |
|
76 NPP = NPP * nsec_per_year |
|
77 |
|
78 ; Redo the units. |
|
79 PREC_ANN@units = "m/yr" |
|
80 RAIN@units = "m/yr" |
|
81 NPP@units = "gC/m^2/yr" |
|
82 TNPP_C@units = "gC/m^2/yr" |
|
83 |
|
84 ;************************************************ |
|
85 pminmax(PREC_ANN,"PREC_ANN") |
|
86 pminmax(TNPP_C,"TNPP_C") |
|
87 pminmax(RAIN,"RAIN") |
|
88 pminmax(NPP,"NPP") |
|
89 |
|
90 RAIN_1D = ndtooned(RAIN) |
|
91 NPP_1D = ndtooned(NPP) |
|
92 PREC_ANN_1D = ndtooned(PREC_ANN) |
|
93 TNPP_C_1D = ndtooned(TNPP_C) |
|
94 |
|
95 ; |
|
96 ; Calculate some "nice" bins for binning the data in equally spaced |
|
97 ; ranges. |
|
98 ; |
|
99 nbins = 15 ; Number of bins to use. |
|
100 |
|
101 nicevals = nice_mnmxintvl(min(RAIN_1D),max(RAIN_1D),nbins,True) |
|
102 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1) |
|
103 range = fspan(nicevals(0),nicevals(1),nvals) |
|
104 ; |
|
105 ; Use this range information to grab all the values in a |
|
106 ; particular range, and then take an average. |
|
107 ; |
|
108 nr = dimsizes(range) |
|
109 nx = nr-1 |
|
110 xvalues = new((/2,nx/),typeof(RAIN_1D)) |
|
111 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. |
|
112 dx = xvalues(0,1) - xvalues(0,0) ; range width |
|
113 dx4 = dx/4 ; 1/4 of the range |
|
114 xvalues(1,:) = xvalues(0,:) - dx/5. |
|
115 yvalues = new((/2,nx/),typeof(RAIN_1D)) |
|
116 mn_yvalues = new((/2,nx/),typeof(RAIN_1D)) |
|
117 mx_yvalues = new((/2,nx/),typeof(RAIN_1D)) |
|
118 |
|
119 do nd=0,1 |
|
120 ; |
|
121 ; See if we are doing model or observational data. |
|
122 ; |
|
123 if(nd.eq.0) then |
|
124 data = PREC_ANN_1D |
|
125 npp_data = TNPP_C_1D |
|
126 else |
|
127 data = RAIN_1D |
|
128 npp_data = NPP_1D |
|
129 end if |
|
130 ; |
|
131 ; Loop through each range and check for values. |
|
132 ; |
|
133 do i=0,nr-2 |
|
134 if (i.ne.(nr-2)) then |
|
135 print("") |
|
136 print("In range ["+range(i)+","+range(i+1)+")") |
|
137 idx = ind((range(i).le.data).and.(data.lt.range(i+1))) |
|
138 else |
|
139 print("") |
|
140 print("In range ["+range(i)+",)") |
|
141 idx = ind(range(i).le.data) |
|
142 end if |
|
143 ; |
|
144 ; Calculate average, and get min and max. |
|
145 ; |
|
146 if(.not.any(ismissing(idx))) then |
|
147 yvalues(nd,i) = avg(npp_data(idx)) |
|
148 mn_yvalues(nd,i) = min(npp_data(idx)) |
|
149 mx_yvalues(nd,i) = max(npp_data(idx)) |
|
150 count = dimsizes(idx) |
|
151 else |
|
152 count = 0 |
|
153 yvalues(nd,i) = yvalues@_FillValue |
|
154 mn_yvalues(nd,i) = yvalues@_FillValue |
|
155 mx_yvalues(nd,i) = yvalues@_FillValue |
|
156 end if |
|
157 ; |
|
158 ; Print out information. |
|
159 ; |
|
160 print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
161 print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
162 |
|
163 ; |
|
164 ; Clean up for next time in loop. |
|
165 ; |
|
166 delete(idx) |
|
167 end do |
|
168 delete(data) |
|
169 delete(npp_data) |
|
170 end do |
|
171 |
|
172 xvalues@long_name = "Mean Annual precipitation (m/year)" |
|
173 yvalues@long_name = "NPP (g C/m2/year)" |
|
174 |
|
175 ; |
|
176 ; Start the graphics. |
|
177 ; |
|
178 wks = gsn_open_wks("png","npp") |
|
179 |
|
180 res = True |
|
181 res@gsnMaximize = False |
|
182 res@gsnDraw = False |
|
183 res@gsnFrame = False |
|
184 res@xyMarkLineMode = "Markers" |
|
185 res@xyMarkerSizeF = 0.014 |
|
186 res@xyMarker = 16 |
|
187 ; res@xyMarkerColors = (/"Gray25","Gray50"/) |
|
188 res@xyMarkerColors = (/"brown","blue"/) |
|
189 res@trYMinF = min(mn_yvalues) - 10. |
|
190 res@trYMaxF = max(mx_yvalues) + 10. |
|
191 ; res@tiMainString = "Observed vs i01.03cn_81site" |
|
192 res@tiMainString = "Observed vs i01.04casa_81site" |
|
193 |
|
194 xy = gsn_csm_xy(wks,xvalues,yvalues,res) |
|
195 |
|
196 max_bar = new((/2,nx/),graphic) |
|
197 min_bar = new((/2,nx/),graphic) |
|
198 max_cap = new((/2,nx/),graphic) |
|
199 min_cap = new((/2,nx/),graphic) |
|
200 |
|
201 lnres = True |
|
202 |
|
203 line_colors = (/"brown","blue"/) |
|
204 do nd=0,1 |
|
205 lnres@gsLineColor = line_colors(nd) |
|
206 do i=0,nx-1 |
|
207 |
|
208 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
209 .not.ismissing(mx_yvalues(nd,i))) then |
|
210 ; |
|
211 ; Attach the vertical bar, both above and below the marker. |
|
212 ; |
|
213 x1 = xvalues(nd,i) |
|
214 y1 = yvalues(nd,i) |
|
215 y2 = mn_yvalues(nd,i) |
|
216 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
217 |
|
218 y2 = mx_yvalues(nd,i) |
|
219 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
220 ; |
|
221 ; Attach the horizontal cap line, both above and below the marker. |
|
222 ; |
|
223 x1 = xvalues(nd,i) - dx4 |
|
224 x2 = xvalues(nd,i) + dx4 |
|
225 y1 = mn_yvalues(nd,i) |
|
226 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
227 |
|
228 y1 = mx_yvalues(nd,i) |
|
229 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
230 end if |
|
231 end do |
|
232 end do |
|
233 |
|
234 draw(xy) |
|
235 frame(wks) |
|
236 |
|
237 end |
|
238 |