Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
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"
8 procedure pminmax(data:numeric,name:string)
10 print ("min/max " + name + " = " + min(data) + "/" + max(data))
11 if(isatt(data,"units")) then
12 print (name + " units = " + data@units)
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")
34 RAIN1 = tofloat(data_file_ob1->LAND_CLASS)
35 NPP1 = data_file_ob2->LAI
36 ;************************************************
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")
44 NPP2 = data_file_model->TLAI
45 ;************************************************
46 ; print min/max and unit
47 ;************************************************
48 pminmax(RAIN1,"RAIN1")
52 RAIN1_1D = ndtooned(RAIN1)
53 NPP1_1D = ndtooned(NPP1)
54 NPP2_1D = ndtooned(NPP2)
56 ; Calculate some "nice" bins for binning the data in equally spaced
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)
66 range = fspan(0,nclassn-1,nclassn)
74 ; Use this range information to grab all the values in a
75 ; particular range, and then take an average.
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))
90 ; See if we are doing model or observational data.
100 ; Loop through each range and check for values.
103 if (i.ne.(nr-2)) then
105 print("In range ["+range(i)+","+range(i+1)+")")
106 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
109 print("In range ["+range(i)+",)")
110 idx = ind(range(i).le.data)
113 ; Calculate average, and get min and max.
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)
122 yvalues(nd,i) = yvalues@_FillValue
123 mn_yvalues(nd,i) = yvalues@_FillValue
124 mx_yvalues(nd,i) = yvalues@_FillValue
127 ; Print out information.
129 print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
130 print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
133 ; Clean up for next time in loop.
142 ; Start the graphics.
150 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
160 bias = sum(abs(vv-uu)/(vv+uu))
161 M = (1.- (bias/nz))*5.