;=======================================================================
; name:
;   xy_zon.ncl
;
; purpose:
;   to plot lon-lat map with zonal mean attached to Y axis
;
; author:
;   m.yoshimori
;=======================================================================

begin

; specify years
  year0 = 1959  ; start year of the dataset
  syear = 1991  ; start year of the analysis
  eyear = 2020  ; end year of the analysis

; specify file names
  fin = "./ERA5/t2m_y1959-y2021.nc"  ; input file name (netCDF)
  fout = "./out/xy_sat_zon"          ; output file name (png/eps etc.)

; get data
  in = addfile(fin,"r")  ; open a file
  tmp = in->var167       ; read the data
  printVarSummary(tmp)   ; check the data

; annual means from "syear" to "eyear"
  y1 = (syear - year0)*12
  y2 = (eyear - year0 + 1)*12 - 1
  yyyymm = yyyymm_time(syear, eyear, "integer")
  tmp2 = month_to_annual_weighted(yyyymm, tmp(y1:y2,:,:), 1)
  delete(tmp)

; long-term average
  var = dim_avg_n_Wrap(tmp2, 0)
  delete(tmp2)

; change unit from (K) to (degC)
  var = var - 273.15

;-----------------
; create base plot
;-----------------

  wks = gsn_open_wks("png",fout)         ; open a png/eps file
  gsn_define_colormap(wks,"temp_19lev")  ; choose a color map

; basic settings
  res = True	              ; plot mods desired
  res@gsnDraw = False         ; draw plot
  res@gsnFrame = False        ; advance frame
  res@cnFillOn = True         ; turn on color fill
  res@cnLinesOn = False       ; turn on contour lines
  res@cnLineLabelsOn = False  ; turn off contour labels

; contour levels
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF = -48.0
  res@cnMaxLevelValF =  48.0
  res@cnLevelSpacingF = 4.0

; adjust aspect ratio
  res@mpShapeMode = "FreeAspect"
  res@vpWidthF = 0.75
  res@vpHeightF = 0.45

; map
  res@mpFillOn = False         ; turn off gray continents
  res@mpCenterLonF = 180       ; set centers the plot at 180
  res@gsnMajorLonSpacing = 60  ; tick mark interval for longitude

; font size
  res@tiMainFontHeightF = 0.022     ; fontsize of the main title
  res@txFontHeightF = 0.02          ; fontsize of the subtitles
  res@tmXBLabelFontHeightF = 0.018  ; fontsize of tickmark labels (x-axis)
  res@tmYLLabelFontHeightF = 0.018  ; fontsize of tickmark labels (x-axis)

; titles
  res@tiMainString = "ERA5 1991-2020 annual mean climatology"  ; title
  res@gsnCenterString = ""                                     ; subtitle
  res@gsnLeftString = "Surface air temperature"                ; upper-left subtitle
  res@gsnRightString = ""                                      ; upper-right subtitle

; colar table
  res@gsnSpreadColors = True  ; spread out color table

; reshape color bar
  res@lbTitleOn          = True
  res@lbOrientation      = "horizontal"
  res@lbBoxEndCapStyle   = "TriangleBothEnds"  ; triangle edges
  res@lbBoxMinorExtentF  = 0.2           ; slim width
  res@lbLabelStride      = 1             ; label interval
  res@pmLabelBarOrthogonalPosF = 0.1   ; shift up&down

; reshape color title
  res@lbTitleString      = "~S~o~N~C"
  res@lbTitleFontHeightF = 0.018
  res@lbTitlePosition    = "Right"
  res@lbTitleDirection   = "Across"

  base_plot = gsn_csm_contour_map(wks,var,res)

;---------------
; create xy plot
;---------------

  x   = dim_avg(var)  ; zonal average
  x!0 = "var"         ; remove warning message
  y   = var&lat       ; create y-axis
  y!0 = "lat"

  xyres                 = True         ; xy plot mods desired
  xyres@gsnDraw         = False        ; don't draw yet
  xyres@gsnFrame        = False        ; don't advance frame yet

  xyres@vpWidthF        = .20          ; set width of second plot
;  xyres@gsnCenterString = "Zonal ave."  ; add title
;  xyres@txFontHeightF   = 0.02          ; change font height

; Bottom axis label
  xyres@tmXBOn = False

; Top axis label
  xyres@tmXTOn = True
  xyres@tmXTLabelsOn = True
  xyres@trXMinF         = -50
  xyres@trXMaxF         = 50
  xyres@trYMinF         = min(var&lat)
  xyres@trYMaxF         = max(var&lat)

  xyres@xyLineThicknesses = 3     ; make lines thicker
  plot = gsn_csm_xy(wks,x,y,xyres)

  x = x*0.
  xyres@xyLineThicknesses = 2     ; make 2nd line thicker
  xyres@xyDashPatterns    = 2     ; make 2nd line dashed
  plot2 = gsn_csm_xy(wks,x,y,xyres)
  overlay(plot,plot2)

;==============================================
; attach plots
;=============================================

  newplot = gsn_attach_plots(base_plot,(/plot/),res,xyres)
  draw(base_plot)
  frame(wks)
 
end