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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

;----------------------------------------------------------------------
; Procedure to attach title to right Y axis
;----------------------------------------------------------------------
procedure add_right_title(wks,plot,title)
local txres, font_height, txid, amres
begin
;---Retrieve font height of left axis string.
  getvalues plot
    "tiXAxisFontHeightF" : font_height
  end getvalues

;---Create a right axis text string to add to plot.
  txres               = True
  txres@txAngleF      = 90.          ; Rotate string clockwise 
  txres@txFontHeightF = font_height  ; Use same font height as left axis

  txid = gsn_create_text(wks, title, txres)


;---Move text string to center/right edge of plot.
  amres                  = True
  amres@amParallelPosF   = 0.68   ; 0.5 is the right edge of the plot, so
                                  ; 0.68 is a little further to the right.
  amres@amOrthogonalPosF = 0.0    ; This is the center of the plot.
  amres@amJust           = "CenterCenter"

  annoid = gsn_add_annotation(plot, txid, amres)  ; Attach string to plot
  return(annoid)
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
 fili  = "mxclim.nc"
 f     = addfile (fili , "r")
 lat   = 30
 t     = f->T(0,:,{lat})    ; Subset T for 1st timestep, all levels, given lat
 printVarSummary(t)

;----------------------------------------------------------------------
; Create "nice height values for labeling the right Y axis later.
;----------------------------------------------------------------------
 hgt   = gsn_geop_hgt(t&lev)    ; This is an unadvertised function.
 hnice = tofloat(ispan(tointeger(floor(hgt(0))), \
         tointeger(ceil(hgt(dimsizes(hgt)-1))),4))
 pnice = ftcurv(hgt,t&lev,hnice)  ; Get pres vals at nice hgt vals.

;----------------------------------------------------------------------
; Start the graphics section
;----------------------------------------------------------------------
 wks   = gsn_open_wks ("x11", "pres_hgt_xy")

;----------------------------------------------------------------------
; First create plot using gsn_csm_pres_hgt with original "u" array
;----------------------------------------------------------------------
 res                 = True

 res@tmYRMode             = "Explicit"
 res@tmYRValues           = pnice
 res@tmYRLabels           = hnice
 res@tmYRMinorOn          = False
 res@tmYRLabelsOn         = True
 res@tmYROn               = True

 res@tmYUseLeft           = False
 res@tmYLMode             = "Explicit"
 res@tmYLMinorOn          = False
 res@tmYLValues           = (/1000, 850, 700, 500, 400, 300, 250,\
                              200,  150, 100,  70,  50,  30,  10/)
 res@tmYLLabels           = "" + res@tmYLValues

 res@trYReverse           = True
 res@trYLog               = True
 res@trYMinF              = min(t&lev)
 res@trYMaxF              = max(t&lev)

 res@tiXAxisString        = ""
 res@tiMainString         = "Temperature (degC) at lat="+lat
 res@tiYAxisString        = "pressure (" + t&lev@units + ")"

 res@gsnDraw              = False
 res@gsnFrame             = False

;---Create plot, add a title on right axis, and draw.
 plot = gsn_csm_xy(wks, t, t&lev, res)
 add_right_title(wks,plot,"height (km)")
 draw(plot)  
 frame(wks)

end