;----------------------------------------------------------------------
; This function is based on the gsn_csm_attach_zonal_means function
;----------------------------------------------------------------------
undef("attach_zonal_plots_to_map_plot")
function attach_zonal_plots_to_map_plot(wks,map_plot,ts1,ts2)
local zres, vph, minlat, maxlat, zplot1, zplot2, anno
begin
;
; Retrieve the lat limits and view port size of the map plot so we plot
; the correct range in Y, and make the plot the correct size.
;
  getvalues map_plot
    "vpHeightF"    : vph
    "mpMinLatF"    : minlat
    "mpMaxLatF"    : maxlat
  end getvalues

  zres                   = True
  zres@gsnDraw           = False
  zres@gsnFrame          = False
  zres@trYMinF           = minlat
  zres@trYMaxF           = maxlat

  zres@trXMinF           = 260   ; I hard-coded these, but you will likely
  zres@trXMaxF           = 310   ; need to change them or better, compute them

  zres@vpWidthF          = 0.12     ; Make the plot thin
  zres@vpHeightF         = vph      ; use same height as map
;  zres@tmYLOn            = False   ; Turns on/off left tickmarks and labels
  zres@tmYROn            = False
  zres@tmXBMaxTicks      = 3
  zres@tmXBMinorPerMajor = 1
  zres@tmXTMinorPerMajor = 1
  zres@tiXAxisOn         = False
  zres@tiYAxisOn         = False

  zres@xyLineThicknessF  = 3.0

  zres@xyLineColor       = "Brown"
  zplot1 = gsn_csm_zonal_means(wks,ts1,zres)
  
  zres@xyLineColor       = "NavyBlue"
  zplot2 = gsn_csm_zonal_means(wks,ts2,zres)

  overlay(zplot1,zplot2)

; Debug
;  draw(zplot1)
;  frame(wks)

;---Add overlaid zonal plots as annotation of map plot.
  anno = NhlAddAnnotation(map_plot,zplot1)
  setvalues anno
    "amZone"          : 2
    "amSide"          : "right"
    "amResizeNotify"  : True
    "amParallelPosF"  : 0.5
    "amOrthogonalPosF": 0.05
  end setvalues

  return(map_plot)
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------

begin
  in  = addfile("ts_Amon_CESM1-CAM5_historical_r1i1p1_185001-200512.nc","r") 
  ts = in->ts
  printVarSummary(ts)

  nt1 = 5      ; pick two arbitary time steps
  nt2 = 1300

  ts1 = ts(nt1,:,:)
  ts2 = ts(nt2,:,:)

  wks  = gsn_open_wks("png","two_zonal_plots")
  
  res                 = True
  
  res@gsnDraw         = False
  res@gsnFrame        = False

  res@vpXF            = 0.01
  res@vpWidthF        = 0.75

  res@cnLineLabelsOn  = False
  res@cnInfoLabelOn   = False

  res@gsnRightString  = ""
  res@gsnLeftString   = ""
  res@tiMainString    = ts@long_name + " (" + ts@units + ")"

  res@trYMinF         = -57
  res@trYMaxF         = 57

  line_plot = gsn_csm_contour(wks,ts1,res)
  delete([/res@trYMinF,res@trYMaxF/])

  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = ispan(210,310,10)

  res@cnFillOn        = True
  res@cnLinesOn       = False

  res@mpProjection    = "Robinson"
  res@mpLimitMode     = "LatLon"
  res@mpMinLatF       = -57
  res@mpMaxLatF       = 57
  res@mpPerimOn       = False

  fill_plot = gsn_csm_contour_map(wks,ts2,res)

  overlay(fill_plot,line_plot)

; Debug
;  draw(fill_plot)
;  frame(wks)

  plot = attach_zonal_plots_to_map_plot(wks,fill_plot,ts1,ts2)

  draw(plot)
  frame(wks)
end