1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/biomass/03.ob_to_T31.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,188 @@
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 +; read in model data
1.15 +;************************************************
1.16 +;fili2 = "b30.061n_1995-2004_ANN_climo_lnd.nc"
1.17 +;model_grid = "T31"
1.18 +;fili2 = "i01.03cn_1545-1569_ANN_climo.nc"
1.19 +;model_grid = "T42"
1.20 + fili2 = "newcn05_ncep_1i_ANN_climo_lnd.nc"
1.21 + model_grid = "1.9"
1.22 +
1.23 + diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.24 + f = addfile (diri2+fili2,"r")
1.25 +
1.26 + lon = f->lon
1.27 + lat = f->lat
1.28 + nlon = dimsizes(lon)
1.29 + nlat = dimsizes(lat)
1.30 +
1.31 +; print (xi)
1.32 +; print (yi)
1.33 +;************************************************
1.34 +; output data
1.35 +;************************************************
1.36 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
1.37 + filo = "amazon_biomass_"+model_grid+".nc"
1.38 + c = addfile(diro+filo,"c")
1.39 +;************************************************
1.40 +; read in observed data
1.41 +;************************************************
1.42 + diri = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
1.43 + fili = "amazon_biomass.nc"
1.44 + g = addfile (diri+fili,"r")
1.45 + bi = g->BIOMASS
1.46 + xi = g->lon
1.47 + yi = g->lat
1.48 +;************************************************
1.49 +; from class to value:
1.50 +; Biomass map has 11 classes: (unit= Mg/ha)
1.51 +; 0: (_FillValue)
1.52 +; 1: 0-25 (12.5)
1.53 +; 2: 25-50 (37.5)
1.54 +; 3: 50-75 (62.5)
1.55 +; 4: 75-100 (87.5)
1.56 +; 5: 100-150 (125)
1.57 +; 6: 150-200 (175)
1.58 +; 7: 200-250 (225)
1.59 +; 8: 250-300 (275)
1.60 +; 9: 300-350 (325)
1.61 +; 10: 350-400 (375)
1.62 +; 11: >400 (425)
1.63 +;--------------------------
1.64 +
1.65 + bi@_FillValue = 1.e+36
1.66 +
1.67 + bi = where(bi.eq.0., bi@_FillValue,bi)
1.68 + bi = where(bi.eq.1., 12.5,bi)
1.69 + bi = where(bi.eq.2., 37.5,bi)
1.70 + bi = where(bi.eq.3., 62.5,bi)
1.71 + bi = where(bi.eq.4., 87.5,bi)
1.72 + bi = where(bi.eq.5., 125.,bi)
1.73 + bi = where(bi.eq.6., 175.,bi)
1.74 + bi = where(bi.eq.7., 225.,bi)
1.75 + bi = where(bi.eq.8., 275.,bi)
1.76 + bi = where(bi.eq.9., 325.,bi)
1.77 + bi = where(bi.eq.10.,375.,bi)
1.78 + bi = where(bi.eq.11.,425.,bi)
1.79 +
1.80 +;print("min/max = " + min(bi) + "/" + max(bi))
1.81 +
1.82 +;************************************************
1.83 +; Observed factor_aboveground = 0.5
1.84 +;************************************************
1.85 +
1.86 + factor_aboveground = 0.5
1.87 +;bi = bi * factor_aboveground
1.88 +
1.89 +
1.90 +;************************************************
1.91 +; find model grids that is inside observed grids
1.92 +;************************************************
1.93 + ind_lonL = min(ind(lon .ge. xi(0)))
1.94 + ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1)))
1.95 +
1.96 + print (ind_lonL)
1.97 +;print (xi(0))
1.98 +;print (lon(ind_lonL))
1.99 +
1.100 + print (ind_lonR)
1.101 +;print (xi(dimsizes(xi)-1))
1.102 +;print (lon(ind_lonR))
1.103 +
1.104 + ind_latS = min(ind(lat .ge. yi(0)))
1.105 + ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1)))
1.106 +
1.107 + print (ind_latS)
1.108 +;print (yi(0))
1.109 +;print (lat(ind_latS))
1.110 +
1.111 + print (ind_latN)
1.112 +;print (yi(dimsizes(yi)-1))
1.113 +;print (lat(ind_latN))
1.114 +
1.115 + nlat_out = ind_latN - ind_latS + 1
1.116 + nlon_out = ind_lonR - ind_lonL + 1
1.117 +
1.118 +print (nlat_out)
1.119 +print (nlon_out)
1.120 +
1.121 + bo = new((/nlat_out,nlon_out/),float)
1.122 + lat_out = new((/nlat_out/),float)
1.123 + lon_out = new((/nlon_out/),float)
1.124 +
1.125 + do jj=0,nlat_out-1
1.126 + j = ind_latS + jj
1.127 + lat_out(jj) = lat(j)
1.128 + if (j.eq.0 .or. j.eq.nlat-1) then
1.129 + if (j.eq.0) then
1.130 + LATS = -90.
1.131 + LATN = lat(j)+0.5*(lat(j+1)-lat(j))
1.132 + end if
1.133 + if (j.eq.nlat-1) then
1.134 + LATS = lat(j)-0.5*(lat(j)-lat(j-1))
1.135 + LATN = 90.
1.136 + end if
1.137 + else
1.138 + LATS = lat(j)-0.5*(lat(j)-lat(j-1))
1.139 + LATN = lat(j)+0.5*(lat(j+1)-lat(j))
1.140 + end if
1.141 + do ii=0,nlon_out-1
1.142 + i = ind_lonL + ii
1.143 + if (jj.eq.0) then
1.144 + lon_out(ii) = lon(i)
1.145 + end if
1.146 + if (i.eq.0 .or. i.eq.nlon-1) then
1.147 + if (i.eq.0) then
1.148 + LONL = 0.
1.149 + LONR = lon(i)+0.5*(lon(i+1)-lon(i))
1.150 + end if
1.151 + if (i.eq.nlon-1) then
1.152 + LONL = lon(i)-0.5*(lon(i)-lon(i-1))
1.153 + LONR = 360.
1.154 + end if
1.155 + else
1.156 + LONL = lon(i)-0.5*(lon(i)-lon(i-1))
1.157 + LONR = lon(i)+0.5*(lon(i+1)-lon(i))
1.158 + end if
1.159 +
1.160 + print (LATS)
1.161 + print (LATN)
1.162 + print (LONL)
1.163 + print (LONR)
1.164 +
1.165 + bo(jj,ii) = avg(bi({LATS:LATN},{LONL:LONR}))
1.166 + end do
1.167 + end do
1.168 +
1.169 + lon_out!0 = "lon"
1.170 + lon_out@long_name = "longitude"
1.171 + lon_out@units = "degrees_east"
1.172 + lon_out&lon = lon_out
1.173 +
1.174 + lat_out!0 = "lat"
1.175 + lat_out@long_name = "latitude"
1.176 + lat_out@units = "degrees_north"
1.177 + lat_out&lat = lat_out
1.178 +
1.179 + bo!0 = "lat"
1.180 + bo!1 = "lon"
1.181 + bo&lat = lat_out
1.182 + bo&lon = lon_out
1.183 + bo@units = bi@units
1.184 + bo@long_name = bi@long_name
1.185 +; bo@_FillValue = bi@_FillValue
1.186 + bo@_FillValue = 1.e+36
1.187 +
1.188 + c->BIOMASS = bo
1.189 + c->lat = lat_out
1.190 + c->lon = lon_out
1.191 +end
1.192 \ No newline at end of file