|
1 ;******************************************************** |
|
2 ; histogram normalized by rain and compute correleration |
|
3 ;******************************************************** |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
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 ; |
|
17 ; Main code. |
|
18 ; |
|
19 begin |
|
20 |
|
21 ;nclass = 18 |
|
22 nclass = 20 |
|
23 |
|
24 ;************************************************ |
|
25 ; read in data: observed |
|
26 ;************************************************ |
|
27 diri1 = "/fis/cgd/cseg/people/jeff/clamp_data/lai/" |
|
28 ;fili1 = "land_class_T42.nc" |
|
29 fili1 = "land_class_T42_new.nc" |
|
30 fili2 = "LAI_2000-2005_mean_T42.nc" |
|
31 data_file_ob1 = addfile(diri1+fili1,"r") |
|
32 data_file_ob2 = addfile(diri1+fili2,"r") |
|
33 |
|
34 RAIN1 = tofloat(data_file_ob1->LAND_CLASS) |
|
35 NPP1 = data_file_ob2->LAI |
|
36 ;************************************************ |
|
37 ; read in data: model |
|
38 ;************************************************ |
|
39 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
40 ;fili3 = "i01.03cn_1545-1569_ANN_climo.nc" |
|
41 fili3 = "i01.04casa_1605-1629_ANN_climo.nc" |
|
42 data_file_model = addfile(diri2+fili3,"r") |
|
43 |
|
44 NPP2 = data_file_model->TLAI |
|
45 ;************************************************ |
|
46 ; print min/max and unit |
|
47 ;************************************************ |
|
48 pminmax(RAIN1,"RAIN1") |
|
49 pminmax(NPP1,"NPP1") |
|
50 pminmax(NPP2,"NPP2") |
|
51 |
|
52 RAIN1_1D = ndtooned(RAIN1) |
|
53 NPP1_1D = ndtooned(NPP1) |
|
54 NPP2_1D = ndtooned(NPP2) |
|
55 ; |
|
56 ; Calculate some "nice" bins for binning the data in equally spaced |
|
57 ; ranges. |
|
58 ; |
|
59 |
|
60 ; nbins = nclass + 1 ; Number of bins to use. |
|
61 ; nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,False) |
|
62 ; nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1) |
|
63 ; range = fspan(nicevals(0),nicevals(1),nvals) |
|
64 |
|
65 nclassn = nclass + 1 |
|
66 range = fspan(0,nclassn-1,nclassn) |
|
67 |
|
68 ; print (nicevals) |
|
69 ; print (nvals) |
|
70 print (range) |
|
71 ; exit |
|
72 |
|
73 ; |
|
74 ; Use this range information to grab all the values in a |
|
75 ; particular range, and then take an average. |
|
76 ; |
|
77 nr = dimsizes(range) |
|
78 nx = nr-1 |
|
79 xvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
80 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. |
|
81 dx = xvalues(0,1) - xvalues(0,0) ; range width |
|
82 dx4 = dx/4 ; 1/4 of the range |
|
83 xvalues(1,:) = xvalues(0,:) - dx/5. |
|
84 yvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
85 mn_yvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
86 mx_yvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
87 |
|
88 do nd=0,1 |
|
89 ; |
|
90 ; See if we are doing model or observational data. |
|
91 ; |
|
92 if(nd.eq.0) then |
|
93 data = RAIN1_1D |
|
94 npp_data = NPP1_1D |
|
95 else |
|
96 data = RAIN1_1D |
|
97 npp_data = NPP2_1D |
|
98 end if |
|
99 ; |
|
100 ; Loop through each range and check for values. |
|
101 ; |
|
102 do i=0,nr-2 |
|
103 if (i.ne.(nr-2)) then |
|
104 print("") |
|
105 print("In range ["+range(i)+","+range(i+1)+")") |
|
106 idx = ind((range(i).le.data).and.(data.lt.range(i+1))) |
|
107 else |
|
108 print("") |
|
109 print("In range ["+range(i)+",)") |
|
110 idx = ind(range(i).le.data) |
|
111 end if |
|
112 ; |
|
113 ; Calculate average, and get min and max. |
|
114 ; |
|
115 if(.not.any(ismissing(idx))) then |
|
116 yvalues(nd,i) = avg(npp_data(idx)) |
|
117 mn_yvalues(nd,i) = min(npp_data(idx)) |
|
118 mx_yvalues(nd,i) = max(npp_data(idx)) |
|
119 count = dimsizes(idx) |
|
120 else |
|
121 count = 0 |
|
122 yvalues(nd,i) = yvalues@_FillValue |
|
123 mn_yvalues(nd,i) = yvalues@_FillValue |
|
124 mx_yvalues(nd,i) = yvalues@_FillValue |
|
125 end if |
|
126 ; |
|
127 ; Print out information. |
|
128 ; |
|
129 print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
130 print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
131 |
|
132 ; |
|
133 ; Clean up for next time in loop. |
|
134 ; |
|
135 delete(idx) |
|
136 end do |
|
137 delete(data) |
|
138 delete(npp_data) |
|
139 end do |
|
140 |
|
141 ; |
|
142 ; Start the graphics. |
|
143 ; |
|
144 |
|
145 u = yvalues(0,:) |
|
146 v = yvalues(1,:) |
|
147 print (u) |
|
148 print (v) |
|
149 |
|
150 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
151 uu = u(good) |
|
152 vv = v(good) |
|
153 nz = dimsizes(uu) |
|
154 print (nz) |
|
155 |
|
156 ccr = esccr(uu,vv,0) |
|
157 print (ccr) |
|
158 |
|
159 ;new eq |
|
160 bias = sum(abs(vv-uu)/(vv+uu)) |
|
161 M = (1.- (bias/nz))*5. |
|
162 print (bias) |
|
163 print (M) |
|
164 |
|
165 end |
|
166 |