;
function yaxpos,x,vec0,vec1,nvec,vb,vt
;             zp   ht
  bracket, x,vec1,nvec,1,k0,k1,ier
  f1 = (x-vec1(k0)) / (vec1(k1)-vec1(k0))
  xvec0 = f1*vec0(k1) + (1.-f1)*vec0(k0)
  ypos = vb + (vt-vb)*(xvec0-vec0(0))/(vec0(nvec-1)-vec0(0))
  return,ypos
end
;-----------------------------------------------------------------------
;
pro bracket, x,xx,nx,inc,n1,n2,ier
;
; Bracket x in xx(nx), returning lower index in n1 and upper index in n2
; If inc > 0 -> array increases from bottom to top,
; If inc <= 0 -> array increases from top to bottom
;
ier = 0
;
; Array increases from bottom to top:
;
if (inc gt 0) then begin
  if (x lt xx(0)) then begin
    n1 = 0
    n2 = 1
    ier = 1
    return
  endif
  if (x gt xx(nx-1)) then begin
    n1 = nx-2
    n2 = nx-1
    ier = 2
    return
  endif
  for i=0,nx-2 do begin
    if (x ge xx(i) and x le xx(i+1)) then begin
      n1 = i
      n2 = i+1
      return
    endif
  endfor
;   
; Array increases from top to bottom:
;
endif else begin
  if (x gt xx(0)) then begin
    n1 = 0
    n2 = 1
    ier = 3
    return
  endif
  if (x lt xx(nx-1)) then begin
    n1 = nx-2
    n2 = nx-1
    ier = 4
    return
  endif
  for i=1,nx-2 do begin
    if (x le xx(i) and x ge xx(i+1)) then begin
      n1 = i
      n2 = i+1
      return
    endif
  endfor
endelse
return
end
;-----------------------------------------------------------------------
pro altyax,info,slice,slicetype,vl,vr,vb,vt,charsize
;this procedure was taken from the tgcmidl model written by Ben Foster
;it was slightly modified to work with the utvert structure
;
zp = *slice.levs 
fields = *info.fields
;
; Get heights from Z field. 
;
IF slicetype eq 'profile' THEN BEGIN
;
;  For profiles only have one Z profile so no need to average
;
  getzprof, info, slice, slicetype, fields, zheights

ENDIF ELSE BEGIN
;
;  Get z coordinate data index
;
  ixz = -1
  for i=0,n_elements(fields)-1 do begin
    if strcompress(fields[i].name,/remove_all) eq 'Z' then ixz = i
    if strcompress(fields[i].name,/remove_all) eq 'Z3' then ixz = i
  endfor
  if ixz eq -1 then begin
    print,'>>> WARNING altyax: Could not find Z field -- extra right-hand',$
      ' axis NOT added.'
    return
  endif
;
; Read Z field if necessary:
; 12/8/05 btf: Currently, altyax is not called for mag fields (see rhyaxis
;              in pltxxx routines). However, ZMAG should be available on
;              histories now, so this routine should check if current field
;              is mag, and if it is, use ZMAG instead of Z.
;
  if not ptr_valid(fields[ixz].data) then begin ; field has not been read
    print,'Reading Z field to get heights for right-hand axes..'
    widget_control,/hourglass
;
; Read Z data from file:
    varget,info,fields[ixz],zdata
    field = fields[ixz]
;
; Process Z field:
    procfield,info,zdata,field,info.z_hist
    fields[ixz] = field
    fmin = min(*fields[ixz].data) & fmax = max(*fields[ixz].data)
    rpt_minmax,info,fields[ixz],fmin,fmax
    *info.fields = fields
  endif
 
  fielddata = *fields[ixz].data
  fieldixz = fields[ixz]
;
;  Need to modify altitude scale for fields at interface levels from WACCM history files.  
;  This is done by using the altitude at midpoints and calculating the values half way 
;  between those points, adding one point near the ground.
;      
  IF (info.ftype EQ 'WACCM' AND slice.field.levstype EQ 'interface') THEN BEGIN 
 
    fieldixz = fields[ixz]
    nLevsIn = fieldixz.nlev
    fileInfo = *info.pfile
    iLevs = fileInfo.ilevs
    nILevs = N_ELEMENTS(iLevs)
 
    inAlts = fielddata

    altsIZPlus1 = inAlts[0:fieldixz.nlon-1,0:fieldixz.nlat-1,1:nLevsIn-1,0:fieldixz.ntime-1]
    altsDiffs = altsIZPlus1 - inAlts[*,*,0:nLevsIn-2,*]
    altsInterface = altsDiffs/2.
    altsILevs = FLTARR(fieldixz.nlon,fieldixz.nlat,nILevs,fieldixz.ntime)
    altsILevs[*,*,0:nILevs-3,*] = inAlts[*,*,0:nLevsIn-2,*] - altsInterface
    altsILevs[*,*,nILevs-2,*] = inAlts[*,*,nLevsIn-1,*] - altsInterface[*,*,nLevsIn-2,*]
    altsILevs[*,*,nILevs-1,*] = inAlts[*,*,nLevsIn-1,*] + altsInterface[*,*,nLevsIn-2,*]
 
    fielddata = altsILevs
    
    print, 'Altitudes converted from midpoint to interface levels '

  ENDIF

  slicedata = *slice.data       ; save
  log10save = slice.log10       ; save
  slice.log10 = 'none'          ; want linear heights
;
;  Pull out 2-D data based on plot type
;
  if slicetype eq 'lon' then begin
    slicename = slice.field.name
    slice.field.name = 'YAXZ'
    deflondata,slice,fielddata   ; get for lon slices (lat,lev)
    slice.field.name = slicename
  endif 
  if slicetype eq 'utvert' then begin
    defutvertdata,slice,fielddata   ; get for utvert slices (ut,lev)
  endif
  if slicetype eq 'lat' then begin 
    deflatdata,slice,fielddata   ; get for lat slices (lon,lev)
  endif
  
  zheights = *slice.data
  *slice.data = slicedata ; restore
;  if slicetype ne 'utvert' then begin
  slice.log10 = log10save ; restore
;  endif

ENDELSE
;
;  For WACCM profiles, the altitude needs to be interpolated 
;
IF info.ftype EQ 'WACCM' AND slicetype EQ 'profile' THEN BEGIN
  
  slicedata = *slice.data
  slice.data = PTR_NEW(zheights)
  interpwdata, slice
  
  zheights = *slice.data
  slice.data = PTR_NEW(slicedata)  

ENDIF
;
; Average over first dimension (lat or lon) except for profiles:
;
IF slicetype eq 'profile' THEN BEGIN

;  izh = WHERE(zp GE slice.levmin-0.01 and zp LE slice.levmax+0.01)
;  nizh = N_ELEMENTS(izh)

  zp = *slice.levs

  IF info.ftype EQ 'WACCM' THEN BEGIN
  
    heights = zheights
    
  ENDIF ELSE BEGIN
;
;  Get the index of the minimum and maximim vertical levels to use in original data  
;  
    levMin = slice.levmin
    iLevMin = IXFIND_NEAREST(zp,levMin)
    levMax = slice.levmax
    iLevMax = IXFIND_NEAREST(zp,levMax)
    
    heights = zheights[iLevMin:iLevMax]
  
  ENDELSE
  
ENDIF ELSE BEGIN

  dims = size(zheights,/dimensions)
  heights = fltarr(dims[1])
  for k=0,dims[1]-1 do begin
    heights[k] = calcmean(zheights[*,k],dims[0],0,1.e-20,$
      float(slice.missing_value),0)
  endfor

ENDELSE

nlev = n_elements(heights)

dht_axis = heights[nlev-1]-heights[0]
;
offtoax  = 0.    ; x offset from existing to new axes
offtonum = .02   ; x offset from existing axis to numeric label on new axis
offtolab = .07   ; x offset from existing axis to new axis string label
offylab  = .009  ; y offset to center numerics next to tic marks
ticlen   = .012  ; major tic len (minor tic len = .5*ticlen)
if charsize eq 0. then chsize = 1.2 else chsize = charsize
;
; Draw new vertical axis offset to right of right axis:
;
xpos = vr+offtoax
plots,[xpos,xpos],[vb,vt],/norm
xyouts,vr+offtolab,.5*(vt+vb),'AVE HEIGHT (KM)',orientation=-90.,/norm,$
  align=.5,charsize=chsize
;
; Draw major tic marks and labels:
;
rndval = 10.
if dht_axis lt 25. then rndval = 5.
if dht_axis lt 20. then rndval = 3.
if dht_axis lt 10. then rndval = 2.
botrnd = rnd(heights[0],rndval,dir=1)       ; round up
toprnd = rnd(heights[nlev-1],rndval,dir=-1) ; round down
nhtlabmin = 6                   ; minimum number of major tics
del = (toprnd-botrnd)/nhtlabmin ; delta
delrnd = rnd(del,rndval)
;if delrnd le 0 then delrnd = 1
if delrnd le 0 then delrnd = toprnd - botrnd
ht = botrnd
while ht le toprnd do begin
  ypos = yaxpos(ht,zp,heights,nlev,vb,vt)
  if ypos le vt and ypos ge vb then begin
    plots,[xpos,xpos+ticlen],[ypos,ypos],/norm
    altlab = strcompress(string(format="(i4)",fix(ht+.5)),/remove_all)
    xyouts,vr+offtonum,ypos-offylab,altlab,/norm,charsize=chsize
  endif
;
; Draw minor tic marks:
;
  if dht_axis ge 10 then begin
    htnext = ht+delrnd                              ; height of next major
    yposnext = yaxpos(htnext,zp,heights,nlev,vb,vt) ; axis position of next major
    nminor = 3 & if yposnext-ypos gt 0.5*(vt-vb) then nminor = 6
    while (htnext-ht) mod nminor ne 0 do nminor=nminor+1
    delminor = (htnext-ht)/nminor
    if delminor le 0 then delminor = 1
    htminor = ht-delrnd
;
; Draw minor tics below bottom major tic:
    while htminor le ht do begin
      yposminor = yaxpos(htminor,zp,heights,nlev,vb,vt)
      if yposminor gt vb then $
        plots,[xpos,xpos+0.5*ticlen],[yposminor,yposminor],/norm
      htminor = htminor+delminor
    end 
    htminor = ht+delminor
;
; Draw minor tics between current and next majors:
    while htminor le htnext do begin
      yposminor = yaxpos(htminor,zp,heights,nlev,vb,vt)
      if yposminor le vt and yposminor ge vb then begin
        plots,[xpos,xpos+0.5*ticlen],[yposminor,yposminor],/norm
      endif 
      if yposnext-ypos gt 0.5*(vt-vb) and htminor eq ht+3.*delminor then begin
        altlab='('+strcompress(string(format="(i4)",fix(ht+.5)),/remove_all)+')'
;       xyouts,vr+offtonum,yposminor-offylab,altlab,/norm,charsize=chsize
      endif
      htminor = htminor+delminor
    end
  endif ; dht_axis ge 10
;
; Increment major tic:
  ht = ht+delrnd
end ; ht le toprnd
end