;*****************************************************
; cru_4.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
;***************************************************
; Specify user desired yyyymm to be plotted
;***************************************************
  yyyymm = (/185101,185107,195001,195007,200001,200007/)
;***************************************************
; Read temperature anomalies
; data are stored as type "short" ... convert to float
;***************************************************
  f      = addfile("crutem2v.nc","r")
  T_anom = short2flt( f->temanom )                   ; contributed.ncl
;***************************************
; For demo purposes: create a companion TIME array use thus to replace the 
; current time coordinate variable, then use this to access specific months
;***************************************  
  time  = f->time            ; read in current array 
  ntime = dimsizes(time)     ; number of months

  TIME = new (ntime,integer)
  delete(TIME@_FillValue)
  TIME@long_name = "TIME"
  TIME@units     = "yyyymm"
  TIME!0         = "TIME"
  
  nmo   = 0
  nyear = 1851
  do nt = 0,ntime-1 
    nmo = nmo+1
    TIME(nt) = nyear*100 + nmo
    if (nmo.eq.12) then
      nmo   = 0
      nyear = nyear+1
    end if
  end do
  TIME&TIME = TIME       ; coordinate variable
  delete (T_anom&time)   ; delete the currently associated coordinate variable
  T_anom&time  = TIME    ; associate the new coordinate variable
;***************************************
; create plot
;***************************************  
  wks = gsn_open_wks("ps","cru")             ; open ps file
  gsn_define_colormap(wks,"gui_default")     ; choose colormap

  plot= new (dimsizes(yyyymm), graphic)      ; create graphical array
  
  res                      = True            ; plot mods desired
  res@gsnDraw              = False           ; don't draw yet 
  res@gsnFrame             = False           ; don't advance frame yet 

  res@cnFillOn             = True            ; color contours
  res@cnLinesOn            = False           ; turn off contour lines
  res@cnFillMode           = "RasterFill"    ; Raster Mode
  res@gsnSpreadColors      = True            ; use full range of color map

  res@lbLabelBarOn         = False           ; turn off individual lb's

  do nt = 0,dimsizes(yyyymm)-1
     res@gsnLeftString   = "Temp Anomaly"                
     res@gsnCenterString = yyyymm(nt)
     plot(nt) = gsn_csm_contour_map_ce(wks,T_anom({yyyymm(nt)},:,:),res)
  end do
;************************************************
; create panel plot
;************************************************
  resP                  = True                   ; modify the panel plot
  resP@gsnPanelLabelBar = True                   ; add common label bar
  gsn_panel(wks,plot,(/3,2/),resP)               ; now draw as one plot

end