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" load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" begin doys = (/"057","089","121","153","185","217"/) ndoy = dimsizes(doys) productid = "MOD15A2" year = 2001 data_dir = "/home/xxi/cdata1/data/modis/"+productid+"/"+year+"/" k = 0 doystr = doys(k) filelist = systemfunc("ls "+data_dir+productid+".A"+year+doystr+"*.hdf") ;print(filelist) nfile = dimsizes(filelist) tiles = str_get_cols(filelist,-28,-23) ;print(tiles) tiles_h = stringtointeger(str_get_cols(tiles,1,2)) tiles_v = stringtointeger(str_get_cols(tiles,4,5)) np = 1200 min_tile_v = min(tiles_v) max_tile_v = max(tiles_v) min_tile_h = min(tiles_h) max_tile_h = max(tiles_h) tvs = ispan(min_tile_v, max_tile_v,1) ths = ispan(min_tile_h, max_tile_h,1) nv = dimsizes(tvs) nh = dimsizes(ths) ;-create coordiantes print("create sinusoidal coordiates") lat2d = new((/nv*np,nh*np/),float) lon2d = new((/nv*np,nh*np/),float) R=6371007.18100 pi=atan(1)*4.0 tile_size=2*pi/36 rad2deg=180.0/pi LineN=conform_dims((/np,np/),ispan(0,np-1,1),0) SampleN=conform_dims((/np,np/),ispan(0,np-1,1),1) do i=0,nv-1 do j=0,nh-1 x = -pi+(ths(j)+(SampleN+0.5)/np)*tile_size y = pi/2-(tvs(i)+(LineN+0.5)/np)*tile_size lat2d(i*np:(i+1)*np-1,j*np:(j+1)*np-1) = y*rad2deg lon2d(i*np:(i+1)*np-1,j*np:(j+1)*np-1) = x/cos(y)*rad2deg end do end do ;- put tiles together lai_all = new((/ndoy,nv*np,nh*np/),float) fpar_all = new((/ndoy,nv*np,nh*np/),float) print("read in data") do k = 0, ndoy-1 lai = new((/nfile,np,np/),float) lai@_FillValue = 999. fpar = new((/nfile,np,np/),float) fpar@_FillValue = 999. do i=0,nfile-1 doystr = doys(k) filelist = systemfunc("ls "+data_dir+productid+".A"+year+doystr+"*.hdf") nfile = dimsizes(filelist) modis_file = addfile(filelist(i), "r") sds_lai = modis_file->Lai_1km ;ubyte, unsigned 8-bit integer sds_fpar = modis_file->Fpar_1km ;ubyte ;sds_qc = modis_file->FparLai_QC ;ubyte tmp1 = tofloat(sds_lai) tmp1@_FillValue = 999. tmp1 = where(tmp1 .ge. 249., tmp1@_FillValue, tmp1) lai(i,:,:) = tmp1 * 0.1 ;scale factor = 0.1 tmp2 = tofloat(sds_fpar) tmp2@_FillValue = 999. tmp2 = where(tmp2 .ge. 100., tmp2@_FillValue, tmp1) fpar(i,:,:) = tmp2 ; unit: % end do do i=0,nfile-1 is = tiles_v(i)-min_tile_v js = tiles_h(i)-min_tile_h lai_all(k,is*np:(is+1)*np-1, js*np:(js+1)*np-1) = lai(i,:,:) fpar_all(k,is*np:(is+1)*np-1, js*np:(js+1)*np-1) = fpar(i,:,:) end do delete(filelist) delete(sds_lai) delete(sds_fpar) delete(lai) delete(fpar) end do lai_all@lat2d=lat2d lai_all@lon2d=lon2d fpar_all@lat2d=lat2d fpar_all@lon2d=lon2d ;- interpolate to wrf domain print("interpolate to wrf domain") wrf_file=addfile("../../thesis/data/geo_usgs_default.nc","r") lat2d_wrf=wrf_file->XLAT_M(0,:,:) lon2d_wrf=wrf_file->XLONG_M(0,:,:) opt = True opt@SrcFileName = "rec_grd.nc" opt@DstFileName = "wrf_grd.nc" opt@WgtFileName = "rec2wrf.nc" opt@SrcGridLat = lat2d ; source grid opt@SrcGridLon = lon2d opt@DstGridLat = lat2d_wrf ; destination grid opt@DstGridLon = lon2d_wrf opt@DstLLCorner = (/min(lat2d_wrf),min(lon2d_wrf)/) opt@DstURCorner = (/max(lat2d_wrf),max(lon2d_wrf)/) opt@SrcRegional = True ; Necessary if grids opt@DstRegional = True ; are regional opt@InterpMethod = "conserve" ; "patch", "conserve" opt@ForceOverwrite = True opt@PrintTimings = True ; Optional. lai_int = ESMF_regrid(lai_all,opt) end