1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/npp/ob/npp_read_avg.ncl	Thu Mar 26 14:31:28 2009 -0400
     1.3 @@ -0,0 +1,88 @@
     1.4 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
     1.5 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
     1.6 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
     1.7 +
     1.8 +begin
     1.9 +  ; Read in binary (2-byte, short int, little endian) MODIS NPP global
    1.10 +  ; annual data at 0.5 degree resolution and create a mean that will
    1.11 +  ; then be reprojected to various computational grids.
    1.12 +
    1.13 +  setfileoption("bin", "ReadByteOrder", "LittleEndian")
    1.14 +  file_prefix = "Npp_0.05deg_"
    1.15 +  file_suffix = ".int16"
    1.16 +  yr_init  = 2000
    1.17 +  ;yr_final = 2005
    1.18 +  yr_final = 2004
    1.19 +
    1.20 +  bin_nlon = 7200
    1.21 +  bin_nlat = 3600
    1.22 +  bin_fill = 32700
    1.23 +  bin_scale_factor = 0.1
    1.24 +
    1.25 +  npp = new((/bin_nlat,bin_nlon/), float)
    1.26 +  npp(:,:) = 0.
    1.27 +
    1.28 +  cnt = 0
    1.29 +  do y = yr_init,yr_final
    1.30 +    ifile = file_prefix + y + file_suffix
    1.31 +    print("Reading " + ifile)
    1.32 +    npp_short = fbindirread(ifile, 0, (/bin_nlat,bin_nlon/), "short")
    1.33 +    npp_short@_FillValue = inttoshort(bin_fill)
    1.34 +    npp = npp + short2flt(npp_short) * bin_scale_factor
    1.35 +    cnt = cnt + 1
    1.36 +  end do
    1.37 +
    1.38 +  delete(npp_short)
    1.39 +
    1.40 +  npp = npp / int2flt(cnt)
    1.41 +
    1.42 +  npp@long_name = "net primary production"
    1.43 +  npp@units     = "gC/m^2/y"
    1.44 +  npp@_FillValue = 1.e20
    1.45 +
    1.46 +  lat = latGlobeFo(bin_nlat, "lat", "latitude", "degrees_north")
    1.47 +  lat = (/ lat(::-1) /)   ; make N->S
    1.48 +  lon = lonGlobeFo(bin_nlon, "lon", "longitude", "degrees_east")
    1.49 +  lon = (/ lon - 180. /)  ; subtract 180 from all values
    1.50 +  lon&lon = lon           ; update coordinates
    1.51 +
    1.52 +  npp!0 = "lat"
    1.53 +  npp!1 = "lon"
    1.54 +  npp&lat = lat
    1.55 +  npp&lon = lon
    1.56 +
    1.57 +  ncfile = addfile("npp_0.05deg_mean_" + yr_init + "-" + yr_final + ".nc", "c")
    1.58 +  ncfile->NPP = npp
    1.59 +
    1.60 +  ;*************************************
    1.61 +  ; create plot
    1.62 +  ;*************************************
    1.63 +
    1.64 +  setvalues NhlGetWorkspaceObjectId()
    1.65 +    "wsMaximumSize" : 199999999
    1.66 +  end setvalues
    1.67 +
    1.68 +  wks = gsn_open_wks("ps", "npp_0.05deg_mean_" + yr_init + "-" + yr_final) ; open a PostScript file
    1.69 +  gsn_define_colormap(wks, "gui_default") ; choose colormap
    1.70 +
    1.71 +  res                    = True         ; Use plot options
    1.72 +  res@cnFillOn           = True         ; Turn on color fill
    1.73 +  res@cnFillMode         = "RasterFill" ; Turn on raster color
    1.74 +  res@lbLabelAutoStride  = True
    1.75 +  res@cnLinesOn          = False        ; Turn off contour lines
    1.76 +  res@gsnSpreadColors    = True         ; Use full colormap
    1.77 +  res@mpFillOn           = False        ; Turn off map fill
    1.78 +
    1.79 +  res@cnLevelSelectionMode = "ManualLevels"    ; Manual contour invtervals
    1.80 +  res@cnMinLevelValF       = 0.                ; Min level
    1.81 +  res@cnMaxLevelValF       = 2200.             ; Max level
    1.82 +  res@cnLevelSpacingF      = 200.              ; interval
    1.83 +  res@tiMainString         = "MODIS mean net primary production (" + yr_init + "-" + yr_final + ")"
    1.84 +
    1.85 +  npp@units     = "gC m~S~-2~N~ y~S~-1~N~"
    1.86 +
    1.87 +  plot = gsn_csm_contour_map_ce(wks,npp,res)
    1.88 +
    1.89 +  delete(wks)
    1.90 +
    1.91 +end