forrest@0: ;******************************************************** forrest@0: ; histogram normalized by rain and compute correleration forrest@0: ;******************************************************** forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" forrest@0: forrest@0: procedure pminmax(data:numeric,name:string) forrest@0: begin forrest@0: print ("min/max " + name + " = " + min(data) + "/" + max(data)) forrest@0: if(isatt(data,"units")) then forrest@0: print (name + " units = " + data@units) forrest@0: end if forrest@0: end forrest@0: forrest@0: ; Main code. forrest@0: begin forrest@0: forrest@0: nclass = 20 forrest@0: forrest@0: plot_type = "ps" forrest@0: plot_type_new = "png" forrest@0: forrest@0: ;************************************************ forrest@0: ; read data: model forrest@0: ;************************************************ forrest@0: co2_i = 283.1878 forrest@0: co2_f = 364.1252 forrest@0: forrest@0: model_grid = "T42" forrest@0: forrest@0: ;model_name_i = "i01.07cn" forrest@0: ;model_name_f = "i01.10cn" forrest@0: forrest@0: model_name_i = "i01.07casa" forrest@0: model_name_f = "i01.10casa" forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film_i = model_name_i + "_1990-2004_ANN_climo.nc" forrest@0: film_f = model_name_f + "_1990-2004_ANN_climo.nc" forrest@0: forrest@0: fm_i = addfile (dirm+film_i,"r") forrest@0: fm_f = addfile (dirm+film_f,"r") forrest@0: forrest@0: npp_i = fm_i->NPP forrest@0: npp_f = fm_f->NPP forrest@0: forrest@0: ;************************************************ forrest@0: ; read data: observed forrest@0: ;************************************************ forrest@0: forrest@0: ob_name = "MODIS MOD 15A2 2000-2005" forrest@0: forrest@0: diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" forrest@0: filo = "land_class_"+model_grid+".nc" forrest@0: forrest@0: fo = addfile(diro+filo,"r") forrest@0: forrest@0: classob = tofloat(fo->LAND_CLASS) forrest@0: forrest@0: class_name = (/"Water Bodies" \ forrest@0: ,"Evergreen Needleleaf Forests" \ forrest@0: ,"Evergreen Broadleaf Forests" \ forrest@0: ,"Deciduous Needleleaf Forest" \ forrest@0: ,"Deciduous Broadleaf Forests" \ forrest@0: ,"Mixed Forests" \ forrest@0: ,"Closed Bushlands" \ forrest@0: ,"Open Bushlands" \ forrest@0: ,"Woody Savannas (S. Hem.)" \ forrest@0: ,"Savannas (S. Hem.)" \ forrest@0: ,"Grasslands" \ forrest@0: ,"Permanent Wetlands" \ forrest@0: ,"Croplands" \ forrest@0: ,"Urban and Built-Up" \ forrest@0: ,"Cropland/Natural Vegetation Mosaic" \ forrest@0: ,"Permanent Snow and Ice" \ forrest@0: ,"Barren or Sparsely Vegetated" \ forrest@0: ,"Unclassified" \ forrest@0: ,"Woody Savannas (N. Hem.)" \ forrest@0: ,"Savannas (N. Hem.)" \ forrest@0: /) forrest@0: forrest@0: ;******************************************************************* forrest@0: ; Calculate "nice" bins for binning the data in equally spaced ranges forrest@0: ;******************************************************************** forrest@0: nclassn = nclass + 1 forrest@0: range = fspan(0,nclassn-1,nclassn) forrest@0: ; print (range) forrest@0: forrest@0: ; Use this range information to grab all the values in a forrest@0: ; particular range, and then take an average. forrest@0: forrest@0: nr = dimsizes(range) forrest@0: nx = nr-1 forrest@0: xvalues = new((/2,nx/),float) forrest@0: xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. forrest@0: dx = xvalues(0,1) - xvalues(0,0) ; range width forrest@0: dx4 = dx/4 ; 1/4 of the range forrest@0: xvalues(1,:) = xvalues(0,:) - dx/5. forrest@0: forrest@0: ; get data forrest@0: forrest@0: DATA11_1D = ndtooned(classob) forrest@0: DATA12_1D = ndtooned(npp_i) forrest@0: DATA22_1D = ndtooned(npp_f) forrest@0: forrest@0: yvalues = new((/2,nx/),float) forrest@0: mn_yvalues = new((/2,nx/),float) forrest@0: mx_yvalues = new((/2,nx/),float) forrest@0: forrest@0: do nd=0,1 forrest@0: forrest@0: ; See if we are doing model or observational data. forrest@0: forrest@0: if(nd.eq.0) then forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA12_1D forrest@0: else forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA22_1D forrest@0: end if forrest@0: forrest@0: ; Loop through each range and check for values. forrest@0: forrest@0: do i=0,nr-2 forrest@0: if (i.ne.(nr-2)) then forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")") forrest@0: idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1))) forrest@0: else forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+",)") forrest@0: idx = ind(data_ob.ge.range(i)) forrest@0: end if forrest@0: forrest@0: ; Calculate average, and get min and max. forrest@0: forrest@0: if(.not.any(ismissing(idx))) then forrest@0: yvalues(nd,i) = avg(data_mod(idx)) forrest@0: mn_yvalues(nd,i) = min(data_mod(idx)) forrest@0: mx_yvalues(nd,i) = max(data_mod(idx)) forrest@0: count = dimsizes(idx) forrest@0: else forrest@0: count = 0 forrest@0: yvalues(nd,i) = yvalues@_FillValue forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue forrest@0: end if forrest@0: forrest@0: ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) forrest@0: forrest@0: ; Clean up for next time in loop. forrest@0: forrest@0: delete(idx) forrest@0: end do forrest@0: delete(data_ob) forrest@0: delete(data_mod) forrest@0: end do forrest@0: forrest@0: ;============================ forrest@0: ;compute beta forrest@0: ;============================ forrest@0: forrest@0: u = yvalues(0,:) forrest@0: v = yvalues(1,:) forrest@0: forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v)) forrest@0: uu = u(good) forrest@0: vv = v(good) forrest@0: ww = class_name(good) forrest@0: forrest@0: n_biome = dimsizes(uu) forrest@0: forrest@0: beta_biome = new((/n_biome/),float) forrest@0: forrest@0: beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i) forrest@0: forrest@0: beta_biome_avg = avg(beta_biome) forrest@0: forrest@0: print("class/beta: " + ww + "/" + beta_biome) forrest@0: print (beta_biome_avg) forrest@0: forrest@0: end forrest@0: