diff -r 000000000000 -r 0c6405ab2ff4 biomass/03.ob_to_T31.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/biomass/03.ob_to_T31.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,188 @@
+; ***********************************************
+; interpolate into model grids (T31)
+; ***********************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+;************************************************
+begin
+
+;************************************************
+; read in model data
+;************************************************
+;fili2  = "b30.061n_1995-2004_ANN_climo_lnd.nc"
+;model_grid = "T31"
+;fili2  = "i01.03cn_1545-1569_ANN_climo.nc"
+;model_grid = "T42"
+ fili2  = "newcn05_ncep_1i_ANN_climo_lnd.nc"
+ model_grid = "1.9"
+
+ diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ f      = addfile (diri2+fili2,"r")
+
+ lon    = f->lon     
+ lat    = f->lat
+ nlon   = dimsizes(lon)
+ nlat   = dimsizes(lat)      
+
+; print (xi)
+; print (yi)
+;************************************************
+; output data
+;************************************************
+  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
+  filo  = "amazon_biomass_"+model_grid+".nc"
+  c = addfile(diro+filo,"c")
+;************************************************
+; read in observed data
+;************************************************
+  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
+  fili  = "amazon_biomass.nc"
+  g     = addfile (diri+fili,"r")
+  bi    = g->BIOMASS   
+  xi    = g->lon 
+  yi    = g->lat
+;************************************************
+; from class to value:
+; Biomass map has 11 classes: (unit= Mg/ha)
+; 0:                (_FillValue)
+; 1:     0-25       (12.5)
+; 2:    25-50       (37.5)
+; 3:    50-75       (62.5)
+; 4:    75-100      (87.5)
+; 5:   100-150      (125)
+; 6:   150-200      (175)
+; 7:   200-250      (225)
+; 8:   250-300      (275)
+; 9:   300-350      (325)
+; 10:  350-400      (375)
+; 11: >400          (425)
+;--------------------------
+
+ bi@_FillValue = 1.e+36
+
+ bi = where(bi.eq.0., bi@_FillValue,bi)
+ bi = where(bi.eq.1., 12.5,bi)
+ bi = where(bi.eq.2., 37.5,bi)
+ bi = where(bi.eq.3., 62.5,bi)
+ bi = where(bi.eq.4., 87.5,bi)
+ bi = where(bi.eq.5., 125.,bi)
+ bi = where(bi.eq.6., 175.,bi)
+ bi = where(bi.eq.7., 225.,bi)
+ bi = where(bi.eq.8., 275.,bi)
+ bi = where(bi.eq.9., 325.,bi)
+ bi = where(bi.eq.10.,375.,bi)
+ bi = where(bi.eq.11.,425.,bi)
+
+;print("min/max = " + min(bi) + "/" + max(bi))
+
+;************************************************
+; Observed factor_aboveground = 0.5
+;************************************************
+
+ factor_aboveground = 0.5
+;bi = bi * factor_aboveground
+
+
+;************************************************
+; find model grids that is inside observed grids
+;************************************************
+ ind_lonL = min(ind(lon .ge. xi(0)))
+ ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1)))
+
+ print (ind_lonL)
+;print (xi(0))
+;print (lon(ind_lonL))
+
+ print (ind_lonR)
+;print (xi(dimsizes(xi)-1))
+;print (lon(ind_lonR))
+
+ ind_latS = min(ind(lat .ge. yi(0)))
+ ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1)))
+
+ print (ind_latS)
+;print (yi(0))
+;print (lat(ind_latS))
+
+ print (ind_latN)
+;print (yi(dimsizes(yi)-1))
+;print (lat(ind_latN))
+
+ nlat_out = ind_latN - ind_latS + 1
+ nlon_out = ind_lonR - ind_lonL + 1
+
+print (nlat_out)
+print (nlon_out)
+ 
+ bo = new((/nlat_out,nlon_out/),float)
+ lat_out = new((/nlat_out/),float)
+ lon_out = new((/nlon_out/),float)
+
+    do jj=0,nlat_out-1
+       j = ind_latS + jj
+       lat_out(jj) = lat(j)
+       if (j.eq.0 .or. j.eq.nlat-1) then
+          if (j.eq.0) then
+             LATS = -90.          
+             LATN = lat(j)+0.5*(lat(j+1)-lat(j))
+          end if
+          if (j.eq.nlat-1) then
+             LATS = lat(j)-0.5*(lat(j)-lat(j-1))
+             LATN = 90.                  
+          end if
+       else
+          LATS = lat(j)-0.5*(lat(j)-lat(j-1))
+          LATN = lat(j)+0.5*(lat(j+1)-lat(j))
+       end if
+       do ii=0,nlon_out-1
+          i = ind_lonL + ii
+          if (jj.eq.0) then
+             lon_out(ii) = lon(i)
+          end if
+          if (i.eq.0 .or. i.eq.nlon-1) then
+             if (i.eq.0) then
+                LONL = 0.          
+                LONR = lon(i)+0.5*(lon(i+1)-lon(i))
+             end if
+             if (i.eq.nlon-1) then
+                LONL = lon(i)-0.5*(lon(i)-lon(i-1))
+                LONR = 360.                 
+             end if
+          else
+             LONL = lon(i)-0.5*(lon(i)-lon(i-1))
+             LONR = lon(i)+0.5*(lon(i+1)-lon(i))
+          end if
+
+ print (LATS)
+ print (LATN)
+ print (LONL)
+ print (LONR)
+
+          bo(jj,ii) = avg(bi({LATS:LATN},{LONL:LONR}))  
+       end do 
+    end do
+
+  lon_out!0          = "lon"
+  lon_out@long_name  = "longitude"
+  lon_out@units      = "degrees_east"
+  lon_out&lon        = lon_out
+
+  lat_out!0          = "lat"
+  lat_out@long_name  = "latitude"
+  lat_out@units      = "degrees_north"
+  lat_out&lat        = lat_out
+
+  bo!0   = "lat"
+  bo!1   = "lon"
+  bo&lat =  lat_out
+  bo&lon =  lon_out
+  bo@units      = bi@units
+  bo@long_name  = bi@long_name
+; bo@_FillValue = bi@_FillValue
+  bo@_FillValue = 1.e+36
+
+  c->BIOMASS  = bo
+  c->lat      = lat_out
+  c->lon      = lon_out
+end
\ No newline at end of file