;---Create dummy time array
function create_time_array(begin_year,end_year)
local units, years, nyears, num_days, ny, nm
begin
  years  = ispan(begin_year,end_year,1)
  nyears = dimsizes(years)
  time   = ispan(0,(nyears*12)-1,1)
  time@units = "months since 2004-01-01"
  return(time)
end

function create_dummy_y_array()
begin
  time = create_time_array(2004,2014)
  ydum = (time%12) + 1
  ydum!0 = "time"
  ydum&time = time
  return(ydum)
end

begin
  y = create_dummy_y_array()

  wks = gsn_open_wks("x11","xy")

  res                   = True

  res@gsnMaximize       = True     ; maximize plot in frame
  res@trYMinF           = min(y)
  res@trXMinF           = min(y&time)
  res@trXMaxF           = max(y&time)

;---Generate XY plot with default X axis labels
  plot = gsn_csm_xy(wks,y&time,y,res)

;---Reverse engineer the time array to get back year, month, day, etc
  dates  = cd_calendar(y&time,0)    ; n x 6 array
  years  = toint(dates(:,0))
  months = toint(dates(:,1))
  days   = toint(dates(:,2))

;---Start constructing the strings for the X axis labels
  month_vals  = (/  3,    6,    9,    12/)     ; months where we want labels
  season_strs = (/"MAM","JJA","SON","DJF"/) + "~C~"   ; ~C~ is a carriage return

;---Get the indexes for the first day of the month for the first month of each season 
  ii  = ind(months.eq.3.or.months.eq.6.or.months.eq.9.or.months.eq.12.and.days.eq.1)
  nii = dimsizes(ii)

  xaxis_labels = new(nii,string)
  do n=0,dimsizes(month_vals)-1
    xaxis_labels = where(months(ii).eq.month_vals(n),season_strs(n)+years(ii),xaxis_labels)
  end do

  res@tmXBMode   = "Explicit"
  res@tmXBValues = y&time(ii(::5))      ; Only label every 4th timestep of the subsetted time array
  res@tmXBLabels = xaxis_labels(::5)
  res@tmXBLabelFontHeightF = 0.012

  plot = gsn_csm_xy(wks,y&time,y,res)
  
end