forrest@0: ; *********************************************** forrest@0: ; interpolate into model grids (T42) forrest@0: ; output: lat: S->N lon: 0->360 forrest@0: ; input : lat: N->S lon: -180->180 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: nclass = 20 forrest@0: forrest@0: ;************************************************ forrest@0: ; output data forrest@0: ;************************************************ forrest@0: diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/" forrest@0: ; filo = "land_class_T42.nc" forrest@0: filo = "land_class_T42_new.nc" forrest@0: c = addfile(diro+filo,"c") forrest@0: forrest@0: ;************************************************ forrest@0: ; read in observed data forrest@0: ;************************************************ forrest@0: diri = "/fis/cgd/cseg/people/jeff/clamp_data/lai/" forrest@0: fili = "land_class_0.25deg.nc" forrest@0: fili = "land_class_0.25deg_new.nc" forrest@0: g = addfile (diri+fili,"r") forrest@0: bi = g->LAND_CLASS forrest@0: xi = g->lon forrest@0: yi = g->lat forrest@0: forrest@0: ;************************************************ forrest@0: ; change into 0-360E, 90S-90N forrest@0: ;************************************************ forrest@0: forrest@0: yi = (/ yi(::-1) /) forrest@0: bi = (/ bi(::-1,:) /) forrest@0: printVarSummary(bi) forrest@0: forrest@0: b2 = bi forrest@0: x2 = xi forrest@0: forrest@0: nx = dimsizes(xi) forrest@0: nx2 = nx/2 forrest@0: forrest@0: do i= 0,nx-1 forrest@0: if (i .lt. nx2) then forrest@0: p = i + nx2 forrest@0: xi(p) = x2(i) + 360. forrest@0: else forrest@0: p = i - nx2 forrest@0: xi(p) = x2(i) forrest@0: end if forrest@0: bi(:,p)= b2(:,i) forrest@0: end do forrest@0: forrest@0: bi&lat = yi forrest@0: bi&lon = xi forrest@0: forrest@0: print (xi) forrest@0: print (yi) forrest@0: ;exit forrest@0: forrest@0: ;************************************************ forrest@0: ; read in model data forrest@0: ;************************************************ forrest@0: diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: fili2 = "i01.03cn_1545-1569_ANN_climo.nc" 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: ; print (xo) forrest@0: ; print (yo) forrest@0: forrest@0: bo = new((/nlat,nlon/),integer) forrest@0: count = new((/nclass/),integer) forrest@0: forrest@0: do j=0,nlat-1 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: forrest@0: do i=0,nlon-1 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: count = 0 forrest@0: forrest@0: do n = 0,nclass-1 forrest@0: count(n) = num (bi({LATS:LATN},{LONL:LONR}).eq.n) forrest@0: end do forrest@0: forrest@0: bo(j,i) = maxind(count) forrest@0: end do forrest@0: end do forrest@0: forrest@0: bo!0 = "lat" forrest@0: bo!1 = "lon" forrest@0: bo&lat = lat forrest@0: bo&lon = lon forrest@0: bo@units = bi@units forrest@0: bo@long_name = bi@long_name forrest@0: bo@_FillValue = bi@_FillValue forrest@0: forrest@0: print (min(bo) + "/" + max(bo)) forrest@0: forrest@0: c->LAND_CLASS = bo forrest@0: end