1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/beta/02.site.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,83 @@
     1.4 +; ****************************************************
     1.5 +; combine  scatter, histogram, global and zonal plots
     1.6 +; *********************************************************
     1.7 +
     1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
     1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    1.11 +; **********************************************************
    1.12 +
    1.13 +; Main code.
    1.14 +begin
    1.15 +
    1.16 +;===================================================
    1.17 +; observed station name and location (lat,lon)
    1.18 +;===================================================
    1.19 +
    1.20 + beta_4_ob = 0.60
    1.21 +
    1.22 + station = (/"DukeFACE" \
    1.23 +            ,"AspenFACE" \
    1.24 +            ,"ORNL-FACE" \
    1.25 +            ,"POP-EUROFACE" \
    1.26 +            /)
    1.27 +
    1.28 + lat_ob = (/ 35.58,  45.40,  35.54, 42.22/)
    1.29 + lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
    1.30 + lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
    1.31 +;print (lon_ob)
    1.32 +
    1.33 + n_sta  = dimsizes(station)
    1.34 +
    1.35 +;====================================================
    1.36 +; model data
    1.37 +;====================================================
    1.38 +
    1.39 + co2_i = 283.1878
    1.40 + co2_f = 364.1252
    1.41 +
    1.42 +;model_name_i = "i01.07cn"
    1.43 +;model_name_f = "i01.10cn"
    1.44 +
    1.45 + model_name_i = "i01.07casa"
    1.46 + model_name_f = "i01.10casa"
    1.47 +
    1.48 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
    1.49 + film_i = model_name_i + "_1990-2004_ANN_climo.nc"
    1.50 + film_f = model_name_f + "_1990-2004_ANN_climo.nc"
    1.51 +
    1.52 + fm_i   = addfile (dirm+film_i,"r")
    1.53 + fm_f   = addfile (dirm+film_f,"r")
    1.54 +  
    1.55 + npp_i  = fm_i->NPP
    1.56 + npp_f  = fm_f->NPP
    1.57 + 
    1.58 + xm       = fm_i->lon  
    1.59 + ym       = fm_i->lat
    1.60 +
    1.61 + npp_i_4  =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
    1.62 +
    1.63 + npp_f_4  =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
    1.64 +
    1.65 +; Units for these variables are:
    1.66 +; npp_i: g C/m^2/s       
    1.67 +
    1.68 + print (npp_i_4)
    1.69 + print (npp_f_4)
    1.70 +;============================
    1.71 +;compute beta
    1.72 +;============================
    1.73 +
    1.74 + beta_4 = new((/n_sta/),float)
    1.75 +
    1.76 + beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
    1.77 +
    1.78 + beta_4_avg = avg(beta_4)
    1.79 +
    1.80 + print (beta_4)
    1.81 + print (beta_4_avg)
    1.82 +
    1.83 + M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 5.
    1.84 + 
    1.85 + print (M_beta)
    1.86 +end