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