forrest@0: ; ****************************************************
forrest@0: ; combine  scatter, histogram, global and zonal plots
forrest@0: ; *********************************************************
forrest@0: 
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ; **********************************************************
forrest@0: 
forrest@0: ; Main code.
forrest@0: begin
forrest@0: 
forrest@0: ;===================================================
forrest@0: ; observed station name and location (lat,lon)
forrest@0: ;===================================================
forrest@0: 
forrest@0:  beta_4_ob = 0.60
forrest@0: 
forrest@0:  station = (/"DukeFACE" \
forrest@0:             ,"AspenFACE" \
forrest@0:             ,"ORNL-FACE" \
forrest@0:             ,"POP-EUROFACE" \
forrest@0:             /)
forrest@0: 
forrest@0:  lat_ob = (/ 35.58,  45.40,  35.54, 42.22/)
forrest@0:  lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
forrest@0:  lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
forrest@0: ;print (lon_ob)
forrest@0: 
forrest@0:  n_sta  = dimsizes(station)
forrest@0: 
forrest@0: ;====================================================
forrest@0: ; model data
forrest@0: ;====================================================
forrest@0: 
forrest@0:  co2_i = 283.1878
forrest@0:  co2_f = 364.1252
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:  xm       = fm_i->lon  
forrest@0:  ym       = fm_i->lat
forrest@0: 
forrest@0:  npp_i_4  =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
forrest@0: 
forrest@0:  npp_f_4  =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
forrest@0: 
forrest@0: ; Units for these variables are:
forrest@0: ; npp_i: g C/m^2/s       
forrest@0: 
forrest@0:  print (npp_i_4)
forrest@0:  print (npp_f_4)
forrest@0: ;============================
forrest@0: ;compute beta
forrest@0: ;============================
forrest@0: 
forrest@0:  beta_4 = new((/n_sta/),float)
forrest@0: 
forrest@0:  beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
forrest@0: 
forrest@0:  beta_4_avg = avg(beta_4)
forrest@0: 
forrest@0:  print (beta_4)
forrest@0:  print (beta_4_avg)
forrest@0: 
forrest@0:  M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 5.
forrest@0:  
forrest@0:  print (M_beta)
forrest@0: end