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