1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lai/41.table_mean.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,166 @@
1.4 +;********************************************************
1.5 +; histogram normalized by rain and compute correleration
1.6 +;********************************************************
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.10 +
1.11 +procedure pminmax(data:numeric,name:string)
1.12 +begin
1.13 + print ("min/max " + name + " = " + min(data) + "/" + max(data))
1.14 + if(isatt(data,"units")) then
1.15 + print (name + " units = " + data@units)
1.16 + end if
1.17 +end
1.18 +
1.19 +;
1.20 +; Main code.
1.21 +;
1.22 +begin
1.23 +
1.24 +;nclass = 18
1.25 + nclass = 20
1.26 +
1.27 +;************************************************
1.28 +; read in data: observed
1.29 +;************************************************
1.30 + diri1 = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
1.31 +;fili1 = "land_class_T42.nc"
1.32 + fili1 = "land_class_T42_new.nc"
1.33 + fili2 = "LAI_2000-2005_mean_T42.nc"
1.34 + data_file_ob1 = addfile(diri1+fili1,"r")
1.35 + data_file_ob2 = addfile(diri1+fili2,"r")
1.36 +
1.37 + RAIN1 = tofloat(data_file_ob1->LAND_CLASS)
1.38 + NPP1 = data_file_ob2->LAI
1.39 +;************************************************
1.40 +; read in data: model
1.41 +;************************************************
1.42 + diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.43 +;fili3 = "i01.03cn_1545-1569_ANN_climo.nc"
1.44 + fili3 = "i01.04casa_1605-1629_ANN_climo.nc"
1.45 + data_file_model = addfile(diri2+fili3,"r")
1.46 +
1.47 + NPP2 = data_file_model->TLAI
1.48 +;************************************************
1.49 +; print min/max and unit
1.50 +;************************************************
1.51 + pminmax(RAIN1,"RAIN1")
1.52 + pminmax(NPP1,"NPP1")
1.53 + pminmax(NPP2,"NPP2")
1.54 +
1.55 + RAIN1_1D = ndtooned(RAIN1)
1.56 + NPP1_1D = ndtooned(NPP1)
1.57 + NPP2_1D = ndtooned(NPP2)
1.58 +;
1.59 +; Calculate some "nice" bins for binning the data in equally spaced
1.60 +; ranges.
1.61 +;
1.62 +
1.63 +; nbins = nclass + 1 ; Number of bins to use.
1.64 +; nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,False)
1.65 +; nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.66 +; range = fspan(nicevals(0),nicevals(1),nvals)
1.67 +
1.68 + nclassn = nclass + 1
1.69 + range = fspan(0,nclassn-1,nclassn)
1.70 +
1.71 +; print (nicevals)
1.72 +; print (nvals)
1.73 + print (range)
1.74 +; exit
1.75 +
1.76 +;
1.77 +; Use this range information to grab all the values in a
1.78 +; particular range, and then take an average.
1.79 +;
1.80 + nr = dimsizes(range)
1.81 + nx = nr-1
1.82 + xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.83 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.84 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.85 + dx4 = dx/4 ; 1/4 of the range
1.86 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.87 + yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.88 + mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.89 + mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.90 +
1.91 + do nd=0,1
1.92 +;
1.93 +; See if we are doing model or observational data.
1.94 +;
1.95 + if(nd.eq.0) then
1.96 + data = RAIN1_1D
1.97 + npp_data = NPP1_1D
1.98 + else
1.99 + data = RAIN1_1D
1.100 + npp_data = NPP2_1D
1.101 + end if
1.102 +;
1.103 +; Loop through each range and check for values.
1.104 +;
1.105 + do i=0,nr-2
1.106 + if (i.ne.(nr-2)) then
1.107 + print("")
1.108 + print("In range ["+range(i)+","+range(i+1)+")")
1.109 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.110 + else
1.111 + print("")
1.112 + print("In range ["+range(i)+",)")
1.113 + idx = ind(range(i).le.data)
1.114 + end if
1.115 +;
1.116 +; Calculate average, and get min and max.
1.117 +;
1.118 + if(.not.any(ismissing(idx))) then
1.119 + yvalues(nd,i) = avg(npp_data(idx))
1.120 + mn_yvalues(nd,i) = min(npp_data(idx))
1.121 + mx_yvalues(nd,i) = max(npp_data(idx))
1.122 + count = dimsizes(idx)
1.123 + else
1.124 + count = 0
1.125 + yvalues(nd,i) = yvalues@_FillValue
1.126 + mn_yvalues(nd,i) = yvalues@_FillValue
1.127 + mx_yvalues(nd,i) = yvalues@_FillValue
1.128 + end if
1.129 +;
1.130 +; Print out information.
1.131 +;
1.132 + print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.133 + print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.134 +
1.135 +;
1.136 +; Clean up for next time in loop.
1.137 +;
1.138 + delete(idx)
1.139 + end do
1.140 + delete(data)
1.141 + delete(npp_data)
1.142 + end do
1.143 +
1.144 +;
1.145 +; Start the graphics.
1.146 +;
1.147 +
1.148 + u = yvalues(0,:)
1.149 + v = yvalues(1,:)
1.150 + print (u)
1.151 + print (v)
1.152 +
1.153 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.154 + uu = u(good)
1.155 + vv = v(good)
1.156 + nz = dimsizes(uu)
1.157 + print (nz)
1.158 +
1.159 + ccr = esccr(uu,vv,0)
1.160 + print (ccr)
1.161 +
1.162 +;new eq
1.163 + bias = sum(abs(vv-uu)/(vv+uu))
1.164 + M = (1.- (bias/nz))*5.
1.165 + print (bias)
1.166 + print (M)
1.167 +
1.168 +end
1.169 +