; ***********************************************

; 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