load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"


begin
;---Read data
  filename  = "WAH_APHRO_TMEAN_1988_2007_MAM_bias.nc"
  a         = addfile(filename,"r")
  rstn_orig = a->rstn

  printVarSummary(rstn_orig)    ; Note the missing value is -99.99
  printMinMax(rstn_orig,0)      ; Note the minimum value is -1.07374e+09

;
; Fix the missing value, since it appears to be incorrect on the file.
;
; There are two ways to do this: set all the current minimum values
; to the missing value defined on the file, or reset the missing value
; to the current minimum value.
;

;---First, get the mininum value and print it.
  rstn_min = min(rstn_orig)   
  print("======================================================================")
  print("Print the rstn minimum without any other arguments to the 'print' call")
  print("  so you can get the full precision.")
  print(rstn_min)
  print("======================================================================")


;---One way to fix the missing value
  rstn_fix = rstn_orig
  rstn_fix@_FillValue = rstn_min

;---Another way to fix the missing value
;  rstn_fix = rstn_orig
;  rstn_fix = where(rstn_orig.eq.rstn_min,rstn_orig@_FillValue,rstn_orig)

  printVarSummary(rstn_fix)
  printMinMax(rstn_fix,0)

;---Just for debug purposes, check how many values we have missing in the fixed data.
  num_missing_orig = num(ismissing(rstn_orig))
  num_missing_fix  = num(ismissing(rstn_fix))
  print("======================================================================")
  print("Number of original values that are equal to the missing value as")
  print("  defined on the file: " + num_missing_orig)
  print("Number of fixed values that are missing: " + num_missing_fix)
  print("======================================================================")

;---Create some graphics
  wks = gsn_open_wks("png","plot_wah")

;---Set some resources
  res                   = True

  res@gsnMaximize       = True     ; maximize plot in frame

  res@cnFillOn          = True     ; turn on contour fill
  res@cnLinesOn         = False    ; turn off contour lines
  res@cnLineLabelsOn    = False    ; turn off line labels
  res@cnFillOpacityF    = 0.5      ; make contour colors half as opaque

  res@tiMainString      = filename
  res@tiMainFontHeightF = 0.02

  res@gsnAddCyclic      = False

;---Zoom in on map area of interest.
  res@mpMinLatF         = min(rstn_orig&lat)-1
  res@mpMaxLatF         = max(rstn_orig&lat)+1
  res@mpMinLonF         = min(rstn_orig&lon)-1
  res@mpMaxLonF         = max(rstn_orig&lon)+1

  res@mpDataBaseVersion     = "MediumRes"      ; better map outlines
  res@mpLandFillColor       = "transparent"    ; default is gray
  res@pmTickMarkDisplayMode = "Always"         ; nicer tickmarks

  plot = gsn_csm_contour_map(wks,rstn_fix(0,:,:),res)

;
; Add markers at each lat/lon location, showing where you have 
; valid data (black) and missing data (red)
;
  mkres                          = True
  mkres@gsMarkerIndex            = 16       ; filled dot
  mkres@gsMarkerSizeF            = 7        ; default is larger
  mkres@gsnCoordsMissingColor    = "red"
  mkres@gsnCoordsNonMissingColor = "black"
  gsn_coordinates(wks,plot,rstn_fix,mkres)

end