clamp: biomass/03.ob_to_T31.ncl@4be95183fbcd

     1 ; ***********************************************

     2 ; interpolate into model grids (T31)

     3 ; ***********************************************

     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

     5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

     7 ;************************************************

     8 begin

    10 ;************************************************

    11 ; read in model data

    12 ;************************************************

    13 ;fili2  = "b30.061n_1995-2004_ANN_climo_lnd.nc"

    14 ;model_grid = "T31"

    15 ;fili2  = "i01.03cn_1545-1569_ANN_climo.nc"

    16 ;model_grid = "T42"

    17  fili2  = "newcn05_ncep_1i_ANN_climo_lnd.nc"

    18  model_grid = "1.9"

    20  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    21  f      = addfile (diri2+fili2,"r")

    23  lon    = f->lon

    24  lat    = f->lat

    25  nlon   = dimsizes(lon)

    26  nlat   = dimsizes(lat)

    28 ; print (xi)

    29 ; print (yi)

    30 ;************************************************

    31 ; output data

    32 ;************************************************

    33   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"

    34   filo  = "amazon_biomass_"+model_grid+".nc"

    35   c = addfile(diro+filo,"c")

    36 ;************************************************

    37 ; read in observed data

    38 ;************************************************

    39   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"

    40   fili  = "amazon_biomass.nc"

    41   g     = addfile (diri+fili,"r")

    42   bi    = g->BIOMASS

    43   xi    = g->lon

    44   yi    = g->lat

    45 ;************************************************

    46 ; from class to value:

    47 ; Biomass map has 11 classes: (unit= Mg/ha)

    48 ; 0:                (_FillValue)

    49 ; 1:     0-25       (12.5)

    50 ; 2:    25-50       (37.5)

    51 ; 3:    50-75       (62.5)

    52 ; 4:    75-100      (87.5)

    53 ; 5:   100-150      (125)

    54 ; 6:   150-200      (175)

    55 ; 7:   200-250      (225)

    56 ; 8:   250-300      (275)

    57 ; 9:   300-350      (325)

    58 ; 10:  350-400      (375)

    59 ; 11: >400          (425)

    60 ;--------------------------

    62  bi@_FillValue = 1.e+36

    64  bi = where(bi.eq.0., bi@_FillValue,bi)

    65  bi = where(bi.eq.1., 12.5,bi)

    66  bi = where(bi.eq.2., 37.5,bi)

    67  bi = where(bi.eq.3., 62.5,bi)

    68  bi = where(bi.eq.4., 87.5,bi)

    69  bi = where(bi.eq.5., 125.,bi)

    70  bi = where(bi.eq.6., 175.,bi)

    71  bi = where(bi.eq.7., 225.,bi)

    72  bi = where(bi.eq.8., 275.,bi)

    73  bi = where(bi.eq.9., 325.,bi)

    74  bi = where(bi.eq.10.,375.,bi)

    75  bi = where(bi.eq.11.,425.,bi)

    77 ;print("min/max = " + min(bi) + "/" + max(bi))

    79 ;************************************************

    80 ; Observed factor_aboveground = 0.5

    81 ;************************************************

    83  factor_aboveground = 0.5

    84 ;bi = bi * factor_aboveground

    87 ;************************************************

    88 ; find model grids that is inside observed grids

    89 ;************************************************

    90  ind_lonL = min(ind(lon .ge. xi(0)))

    91  ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1)))

    93  print (ind_lonL)

    94 ;print (xi(0))

    95 ;print (lon(ind_lonL))

    97  print (ind_lonR)

    98 ;print (xi(dimsizes(xi)-1))

    99 ;print (lon(ind_lonR))

   101  ind_latS = min(ind(lat .ge. yi(0)))

   102  ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1)))

   104  print (ind_latS)

   105 ;print (yi(0))

   106 ;print (lat(ind_latS))

   108  print (ind_latN)

   109 ;print (yi(dimsizes(yi)-1))

   110 ;print (lat(ind_latN))

   112  nlat_out = ind_latN - ind_latS + 1

   113  nlon_out = ind_lonR - ind_lonL + 1

   115 print (nlat_out)

   116 print (nlon_out)

   118  bo = new((/nlat_out,nlon_out/),float)

   119  lat_out = new((/nlat_out/),float)

   120  lon_out = new((/nlon_out/),float)

   122     do jj=0,nlat_out-1

   123        j = ind_latS + jj

   124        lat_out(jj) = lat(j)

   125        if (j.eq.0 .or. j.eq.nlat-1) then

   126           if (j.eq.0) then

   127              LATS = -90.

   128              LATN = lat(j)+0.5*(lat(j+1)-lat(j))

   129           end if

   130           if (j.eq.nlat-1) then

   131              LATS = lat(j)-0.5*(lat(j)-lat(j-1))

   132              LATN = 90.

   133           end if

   134        else

   135           LATS = lat(j)-0.5*(lat(j)-lat(j-1))

   136           LATN = lat(j)+0.5*(lat(j+1)-lat(j))

   137        end if

   138        do ii=0,nlon_out-1

   139           i = ind_lonL + ii

   140           if (jj.eq.0) then

   141              lon_out(ii) = lon(i)

   142           end if

   143           if (i.eq.0 .or. i.eq.nlon-1) then

   144              if (i.eq.0) then

   145                 LONL = 0.

   146                 LONR = lon(i)+0.5*(lon(i+1)-lon(i))

   147              end if

   148              if (i.eq.nlon-1) then

   149                 LONL = lon(i)-0.5*(lon(i)-lon(i-1))

   150                 LONR = 360.

   151              end if

   152           else

   153              LONL = lon(i)-0.5*(lon(i)-lon(i-1))

   154              LONR = lon(i)+0.5*(lon(i+1)-lon(i))

   155           end if

   157  print (LATS)

   158  print (LATN)

   159  print (LONL)

   160  print (LONR)

   162           bo(jj,ii) = avg(bi({LATS:LATN},{LONL:LONR}))

   163        end do

   164     end do

   166   lon_out!0          = "lon"

   167   lon_out@long_name  = "longitude"

   168   lon_out@units      = "degrees_east"

   169   lon_out&lon        = lon_out

   171   lat_out!0          = "lat"

   172   lat_out@long_name  = "latitude"

   173   lat_out@units      = "degrees_north"

   174   lat_out&lat        = lat_out

   176   bo!0   = "lat"

   177   bo!1   = "lon"

   178   bo&lat =  lat_out

   179   bo&lon =  lon_out

   180   bo@units      = bi@units

   181   bo@long_name  = bi@long_name

   182 ; bo@_FillValue = bi@_FillValue

   183   bo@_FillValue = 1.e+36

   185   c->BIOMASS  = bo

   186   c->lat      = lat_out

   187   c->lon      = lon_out

   188 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--