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"

 begin

   ; Read in binary (2-byte, short int, little endian) MODIS NPP global

   ; annual data at 0.5 degree resolution and create a mean that will

   ; then be reprojected to various computational grids.

   setfileoption("bin", "ReadByteOrder", "LittleEndian")

   file_prefix = "Npp_0.05deg_"

   file_suffix = ".int16"

   yr_init  = 2000

   ;yr_final = 2005

   yr_final = 2004

   bin_nlon = 7200

   bin_nlat = 3600

   bin_fill = 32700

   bin_scale_factor = 0.1

   npp = new((/bin_nlat,bin_nlon/), float)

   npp(:,:) = 0.

   cnt = 0

   do y = yr_init,yr_final

     ifile = file_prefix + y + file_suffix

     print("Reading " + ifile)

     npp_short = fbindirread(ifile, 0, (/bin_nlat,bin_nlon/), "short")

     npp_short@_FillValue = inttoshort(bin_fill)

     npp = npp + short2flt(npp_short) * bin_scale_factor

     cnt = cnt + 1

   end do

   delete(npp_short)

   npp = npp / int2flt(cnt)

   npp@long_name = "net primary production"

   npp@units     = "gC/m^2/y"

   npp@_FillValue = 1.e20

   lat = latGlobeFo(bin_nlat, "lat", "latitude", "degrees_north")

   lat = (/ lat(::-1) /)   ; make N->S

   lon = lonGlobeFo(bin_nlon, "lon", "longitude", "degrees_east")

   lon = (/ lon - 180. /)  ; subtract 180 from all values

   lon&lon = lon           ; update coordinates

   npp!0 = "lat"

   npp!1 = "lon"

   npp&lat = lat

   npp&lon = lon

   ncfile = addfile("npp_0.05deg_mean_" + yr_init + "-" + yr_final + ".nc", "c")

   ncfile->NPP = npp

   ;*************************************

   ; create plot

   ;*************************************

   setvalues NhlGetWorkspaceObjectId()

     "wsMaximumSize" : 199999999

   end setvalues

   wks = gsn_open_wks("ps", "npp_0.05deg_mean_" + yr_init + "-" + yr_final) ; open a PostScript file

   gsn_define_colormap(wks, "gui_default") ; choose colormap

   res                    = True         ; Use plot options

   res@cnFillOn           = True         ; Turn on color fill

   res@cnFillMode         = "RasterFill" ; Turn on raster color

   res@lbLabelAutoStride  = True

   res@cnLinesOn          = False        ; Turn off contour lines

   res@gsnSpreadColors    = True         ; Use full colormap

   res@mpFillOn           = False        ; Turn off map fill

   res@cnLevelSelectionMode = "ManualLevels"    ; Manual contour invtervals

   res@cnMinLevelValF       = 0.                ; Min level

   res@cnMaxLevelValF       = 2200.             ; Max level

   res@cnLevelSpacingF      = 200.              ; interval

   res@tiMainString         = "MODIS mean net primary production (" + yr_init + "-" + yr_final + ")"

   npp@units     = "gC m~S~-2~N~ y~S~-1~N~"

   plot = gsn_csm_contour_map_ce(wks,npp,res)

   delete(wks)

 end