;***********************************************
; coads_3.ncl
;***********************************************
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
;*****************************************************
; get data
;*****************************************************
  fa     = addfile("air.mean.nc","r")
  fb     = addfile("air.nobs.nc","r")
 
  air  = short2flt( fa->air)                           ; air temp        
  nobs = short2flt( fb->air)                           ; # obs for air
; filter out geographic locations that have < nMin observations
  nMin   = 2                                           ; min # obs
  air    = mask(air,nobs.ge.nMin,True)
;************************************************
; Compute the climatology using a function in contributed.ncl
;************************************************
  yrStrt = 1950
  yrLast = 1979
  moStrt = (yrStrt-1800)*12              ; start subscript
  moLast = (yrLast-1800)*12 + 11         ; last  subscript
  clm = clmMonTLL (air(moStrt:moLast,:,:)) ; monthly climatology
;************************************************
; Compute the anomalies from the climatology just computed
;************************************************
  xAnom = calcMonAnomTLL (air,clm)
;************************************************
; Graphics
;*************************************************
  wks = gsn_open_wks("ps","coads")                ; open ps file  
  gsn_define_colormap(wks,"BlWhRe")               ; choose colormap
  i = NhlNewColor(wks,0.8,0.8,0.8)                ; add gray to colormap 

  res                      = True                 ; plot options desired
  res@cnFillOn             = True                 ; turn on color fill
  res@cnLinesOn            = False                ; turn off contour lines
  res@cnLineLabelsOn       = False                ; turn off line labels

  res@gsnSpreadColors      = True                 ; use full color map
  res@gsnSpreadColorStart  = 32                   ; skip dark blues
  res@gsnSpreadColorEnd    = 72                   ; skip dark reds

  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       =  -3.               ; set min contour level
  res@cnMaxLevelValF       =   3.               ; set max contour level
  res@cnLevelSpacingF      =   0.5              ; set contour spacing

  res@gsnMaximize          = True                 ; blow up plot
      
  nt                       = 2200   
  res@gsnRightString       = "time="+xAnom&time(nt)  
  plot  = gsn_csm_contour_map_ce(wks,xAnom(nt,:,:),res) 
end