1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/42.bias_zonal.nc Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,42 @@
1.4 +; ***********************************************
1.5 +; xy_4.ncl
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 +begin
1.12 +;************************************************
1.13 +; read in observed data
1.14 +;************************************************
1.15 + f = addfile ("Npp_T42_mean.nc","r")
1.16 + y = f->NPP
1.17 +;************************************************
1.18 +; read in model data
1.19 +;************************************************
1.20 + g = addfile ("i01.03cn_1545-1569_ANN_climo.nc","r")
1.21 +;g = addfile ("i01.04casa_1605-1629_ANN_climo.nc","r")
1.22 + x = g->NPP
1.23 +
1.24 + delta = 0.00000000001
1.25 + x0 = x(0,:,:)
1.26 + y = where(ismissing(y).and.(ismissing(x0).or.(x0.lt.delta)),0.,y)
1.27 +
1.28 + p = zonalAve(y)
1.29 +
1.30 + sec_to_year = 86400.*365.
1.31 + x0 = x0 * sec_to_year
1.32 + q = zonalAve(x0)
1.33 +
1.34 + good = ind(p .ne. 0.)
1.35 + u = p(good)
1.36 + v = q(good)
1.37 +
1.38 + ccr = esccr(u,v,0)
1.39 + print (ccr)
1.40 + bias = sum(((v-u)/u)^2)
1.41 + print (bias)
1.42 + M = 1.- sqrt(bias/dimsizes(u))
1.43 + print (M)
1.44 +
1.45 +end