1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/biomass/03.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,115 @@
1.4 +; ***********************************************
1.5 +; interpolate into model grids (T31)
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 +
1.14 +;************************************************
1.15 +; read in observed data
1.16 +;************************************************
1.17 + diri = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
1.18 + fili = "amazon_biomass.nc"
1.19 + g = addfile (diri+fili,"r")
1.20 + bi = g->BIOMASS
1.21 + xi = g->lon
1.22 + yi = g->lat
1.23 +
1.24 +
1.25 +;************************************************
1.26 +; read in model data
1.27 +;************************************************
1.28 + diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.29 + fili2 = "b30.061n_1995-2004_ANN_climo_lnd.nc"
1.30 + f = addfile (diri2+fili2,"r")
1.31 +
1.32 + lon = f->lon
1.33 + lat = f->lat
1.34 + nlon = dimsizes(lon)
1.35 + nlat = dimsizes(lat)
1.36 +
1.37 +; print (xi)
1.38 +; print (yi)
1.39 +
1.40 + ind_lonL = min(ind(lon .ge. xi(0)))
1.41 + ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1)))
1.42 +
1.43 +print (ind_lonL)
1.44 +print (xi(0))
1.45 +print (lon(ind_lonL))
1.46 +
1.47 +print (ind_lonR)
1.48 +print (xi(dimsizes(xi)-1))
1.49 +print (lon(ind_lonR))
1.50 +
1.51 + ind_latS = min(ind(lat .ge. yi(0)))
1.52 + ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1)))
1.53 +
1.54 +print (ind_latS)
1.55 +print (yi(0))
1.56 +print (lat(ind_latS))
1.57 +
1.58 +print (ind_latN)
1.59 +print (yi(dimsizes(yi)-1))
1.60 +print (lat(ind_latN))
1.61 +exit
1.62 + bo = new((/nlat,nlon/),float)
1.63 + do j=0,nlat-1
1.64 + if (j.eq.0 .or. j.eq.nlat-1) then
1.65 + if (j.eq.0) then
1.66 + LATS = -90.
1.67 + LATN = lat(j)+0.5*(lat(j+1)-lat(j))
1.68 + end if
1.69 + if (j.eq.nlat-1) then
1.70 + LATS = lat(j)-0.5*(lat(j)-lat(j-1))
1.71 + LATN = 90.
1.72 + end if
1.73 + else
1.74 + LATS = lat(j)-0.5*(lat(j)-lat(j-1))
1.75 + LATN = lat(j)+0.5*(lat(j+1)-lat(j))
1.76 + end if
1.77 +
1.78 +; CLAT = clat({LATS:LATN}) ; do once for *slight* efficiency
1.79 +; TEMP = bi({LATS:LATN},:) ; 2D [lat,lon]
1.80 +
1.81 + do i=0,nlon-1
1.82 + if (i.eq.0 .or. i.eq.nlon-1) then
1.83 + if (i.eq.0) then
1.84 + LONL = 0.
1.85 + LONR = lon(i)+0.5*(lon(i+1)-lon(i))
1.86 + end if
1.87 + if (i.eq.nlon-1) then
1.88 + LONL = lon(i)-0.5*(lon(i)-lon(i-1))
1.89 + LONR = 360.
1.90 + end if
1.91 + else
1.92 + LONL = lon(i)-0.5*(lon(i)-lon(i-1))
1.93 + LONR = lon(i)+0.5*(lon(i+1)-lon(i))
1.94 + end if
1.95 +
1.96 +;print (LATS)
1.97 +;print (LATN)
1.98 +;print (LONL)
1.99 +;print (LONR)
1.100 + bo(j,i) = avg(bi({LATS:LATN},{LONL:LONR}))
1.101 +; bo(j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
1.102 + end do
1.103 +
1.104 +; delete(CLAT)
1.105 +; delete(TEMP)
1.106 + end do
1.107 +
1.108 + bo!0 = "lat"
1.109 + bo!1 = "lon"
1.110 + bo&lat = lat
1.111 + bo&lon = lon
1.112 + bo@units = bi@units
1.113 + bo@long_name = bi@long_name
1.114 +; bo@_FillValue = bi@_FillValue
1.115 + bo@_FillValue = 1.e+36
1.116 +
1.117 + c->NPP = bo
1.118 +end
1.119 \ No newline at end of file