1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/beta/03.biome.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,185 @@
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.test"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
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 +; Main code.
1.20 +begin
1.21 +
1.22 + nclass = 20
1.23 +
1.24 + plot_type = "ps"
1.25 + plot_type_new = "png"
1.26 +
1.27 +;************************************************
1.28 +; read data: model
1.29 +;************************************************
1.30 + co2_i = 283.1878
1.31 + co2_f = 364.1252
1.32 +
1.33 + model_grid = "T42"
1.34 +
1.35 +;model_name_i = "i01.07cn"
1.36 +;model_name_f = "i01.10cn"
1.37 +
1.38 + model_name_i = "i01.07casa"
1.39 + model_name_f = "i01.10casa"
1.40 +
1.41 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.42 + film_i = model_name_i + "_1990-2004_ANN_climo.nc"
1.43 + film_f = model_name_f + "_1990-2004_ANN_climo.nc"
1.44 +
1.45 + fm_i = addfile (dirm+film_i,"r")
1.46 + fm_f = addfile (dirm+film_f,"r")
1.47 +
1.48 + npp_i = fm_i->NPP
1.49 + npp_f = fm_f->NPP
1.50 +
1.51 +;************************************************
1.52 +; read data: observed
1.53 +;************************************************
1.54 +
1.55 + ob_name = "MODIS MOD 15A2 2000-2005"
1.56 +
1.57 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
1.58 + filo = "land_class_"+model_grid+".nc"
1.59 +
1.60 + fo = addfile(diro+filo,"r")
1.61 +
1.62 + classob = tofloat(fo->LAND_CLASS)
1.63 +
1.64 + class_name = (/"Water Bodies" \
1.65 + ,"Evergreen Needleleaf Forests" \
1.66 + ,"Evergreen Broadleaf Forests" \
1.67 + ,"Deciduous Needleleaf Forest" \
1.68 + ,"Deciduous Broadleaf Forests" \
1.69 + ,"Mixed Forests" \
1.70 + ,"Closed Bushlands" \
1.71 + ,"Open Bushlands" \
1.72 + ,"Woody Savannas (S. Hem.)" \
1.73 + ,"Savannas (S. Hem.)" \
1.74 + ,"Grasslands" \
1.75 + ,"Permanent Wetlands" \
1.76 + ,"Croplands" \
1.77 + ,"Urban and Built-Up" \
1.78 + ,"Cropland/Natural Vegetation Mosaic" \
1.79 + ,"Permanent Snow and Ice" \
1.80 + ,"Barren or Sparsely Vegetated" \
1.81 + ,"Unclassified" \
1.82 + ,"Woody Savannas (N. Hem.)" \
1.83 + ,"Savannas (N. Hem.)" \
1.84 + /)
1.85 +
1.86 +;*******************************************************************
1.87 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.88 +;********************************************************************
1.89 + nclassn = nclass + 1
1.90 + range = fspan(0,nclassn-1,nclassn)
1.91 +; print (range)
1.92 +
1.93 +; Use this range information to grab all the values in a
1.94 +; particular range, and then take an average.
1.95 +
1.96 + nr = dimsizes(range)
1.97 + nx = nr-1
1.98 + xvalues = new((/2,nx/),float)
1.99 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.100 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.101 + dx4 = dx/4 ; 1/4 of the range
1.102 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.103 +
1.104 +; get data
1.105 +
1.106 + DATA11_1D = ndtooned(classob)
1.107 + DATA12_1D = ndtooned(npp_i)
1.108 + DATA22_1D = ndtooned(npp_f)
1.109 +
1.110 + yvalues = new((/2,nx/),float)
1.111 + mn_yvalues = new((/2,nx/),float)
1.112 + mx_yvalues = new((/2,nx/),float)
1.113 +
1.114 + do nd=0,1
1.115 +
1.116 +; See if we are doing model or observational data.
1.117 +
1.118 + if(nd.eq.0) then
1.119 + data_ob = DATA11_1D
1.120 + data_mod = DATA12_1D
1.121 + else
1.122 + data_ob = DATA11_1D
1.123 + data_mod = DATA22_1D
1.124 + end if
1.125 +
1.126 +; Loop through each range and check for values.
1.127 +
1.128 + do i=0,nr-2
1.129 + if (i.ne.(nr-2)) then
1.130 +; print("")
1.131 +; print("In range ["+range(i)+","+range(i+1)+")")
1.132 + idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))
1.133 + else
1.134 +; print("")
1.135 +; print("In range ["+range(i)+",)")
1.136 + idx = ind(data_ob.ge.range(i))
1.137 + end if
1.138 +
1.139 +; Calculate average, and get min and max.
1.140 +
1.141 + if(.not.any(ismissing(idx))) then
1.142 + yvalues(nd,i) = avg(data_mod(idx))
1.143 + mn_yvalues(nd,i) = min(data_mod(idx))
1.144 + mx_yvalues(nd,i) = max(data_mod(idx))
1.145 + count = dimsizes(idx)
1.146 + else
1.147 + count = 0
1.148 + yvalues(nd,i) = yvalues@_FillValue
1.149 + mn_yvalues(nd,i) = yvalues@_FillValue
1.150 + mx_yvalues(nd,i) = yvalues@_FillValue
1.151 + end if
1.152 +
1.153 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.154 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.155 +
1.156 +; Clean up for next time in loop.
1.157 +
1.158 + delete(idx)
1.159 + end do
1.160 + delete(data_ob)
1.161 + delete(data_mod)
1.162 + end do
1.163 +
1.164 +;============================
1.165 +;compute beta
1.166 +;============================
1.167 +
1.168 + u = yvalues(0,:)
1.169 + v = yvalues(1,:)
1.170 +
1.171 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.172 + uu = u(good)
1.173 + vv = v(good)
1.174 + ww = class_name(good)
1.175 +
1.176 + n_biome = dimsizes(uu)
1.177 +
1.178 + beta_biome = new((/n_biome/),float)
1.179 +
1.180 + beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
1.181 +
1.182 + beta_biome_avg = avg(beta_biome)
1.183 +
1.184 + print("class/beta: " + ww + "/" + beta_biome)
1.185 + print (beta_biome_avg)
1.186 +
1.187 +end
1.188 +