npp/ob/npp_read_avg.ncl
changeset 2 e7ba9bcc3020
equal deleted inserted replaced
-1:000000000000 0:ae94a6ce7583
       
     1 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
     2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
     3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
     4 
       
     5 begin
       
     6   ; Read in binary (2-byte, short int, little endian) MODIS NPP global
       
     7   ; annual data at 0.5 degree resolution and create a mean that will
       
     8   ; then be reprojected to various computational grids.
       
     9 
       
    10   setfileoption("bin", "ReadByteOrder", "LittleEndian")
       
    11   file_prefix = "Npp_0.05deg_"
       
    12   file_suffix = ".int16"
       
    13   yr_init  = 2000
       
    14   ;yr_final = 2005
       
    15   yr_final = 2004
       
    16 
       
    17   bin_nlon = 7200
       
    18   bin_nlat = 3600
       
    19   bin_fill = 32700
       
    20   bin_scale_factor = 0.1
       
    21 
       
    22   npp = new((/bin_nlat,bin_nlon/), float)
       
    23   npp(:,:) = 0.
       
    24 
       
    25   cnt = 0
       
    26   do y = yr_init,yr_final
       
    27     ifile = file_prefix + y + file_suffix
       
    28     print("Reading " + ifile)
       
    29     npp_short = fbindirread(ifile, 0, (/bin_nlat,bin_nlon/), "short")
       
    30     npp_short@_FillValue = inttoshort(bin_fill)
       
    31     npp = npp + short2flt(npp_short) * bin_scale_factor
       
    32     cnt = cnt + 1
       
    33   end do
       
    34 
       
    35   delete(npp_short)
       
    36 
       
    37   npp = npp / int2flt(cnt)
       
    38 
       
    39   npp@long_name = "net primary production"
       
    40   npp@units     = "gC/m^2/y"
       
    41   npp@_FillValue = 1.e20
       
    42 
       
    43   lat = latGlobeFo(bin_nlat, "lat", "latitude", "degrees_north")
       
    44   lat = (/ lat(::-1) /)   ; make N->S
       
    45   lon = lonGlobeFo(bin_nlon, "lon", "longitude", "degrees_east")
       
    46   lon = (/ lon - 180. /)  ; subtract 180 from all values
       
    47   lon&lon = lon           ; update coordinates
       
    48 
       
    49   npp!0 = "lat"
       
    50   npp!1 = "lon"
       
    51   npp&lat = lat
       
    52   npp&lon = lon
       
    53 
       
    54   ncfile = addfile("npp_0.05deg_mean_" + yr_init + "-" + yr_final + ".nc", "c")
       
    55   ncfile->NPP = npp
       
    56 
       
    57   ;*************************************
       
    58   ; create plot
       
    59   ;*************************************
       
    60 
       
    61   setvalues NhlGetWorkspaceObjectId()
       
    62     "wsMaximumSize" : 199999999
       
    63   end setvalues
       
    64 
       
    65   wks = gsn_open_wks("ps", "npp_0.05deg_mean_" + yr_init + "-" + yr_final) ; open a PostScript file
       
    66   gsn_define_colormap(wks, "gui_default") ; choose colormap
       
    67 
       
    68   res                    = True         ; Use plot options
       
    69   res@cnFillOn           = True         ; Turn on color fill
       
    70   res@cnFillMode         = "RasterFill" ; Turn on raster color
       
    71   res@lbLabelAutoStride  = True
       
    72   res@cnLinesOn          = False        ; Turn off contour lines
       
    73   res@gsnSpreadColors    = True         ; Use full colormap
       
    74   res@mpFillOn           = False        ; Turn off map fill
       
    75 
       
    76   res@cnLevelSelectionMode = "ManualLevels"    ; Manual contour invtervals
       
    77   res@cnMinLevelValF       = 0.                ; Min level
       
    78   res@cnMaxLevelValF       = 2200.             ; Max level
       
    79   res@cnLevelSpacingF      = 200.              ; interval
       
    80   res@tiMainString         = "MODIS mean net primary production (" + yr_init + "-" + yr_final + ")"
       
    81 
       
    82   npp@units     = "gC m~S~-2~N~ y~S~-1~N~"
       
    83 
       
    84   plot = gsn_csm_contour_map_ce(wks,npp,res)
       
    85 
       
    86   delete(wks)
       
    87 
       
    88 end